Abstract
We consider the numerical integration of a size-structured cell population model. We propose a new second-order numerical method to attain its solution. The scheme is analyzed and the optimal rate of convergence is derived. We show experimentally the predicted accuracy of the scheme.
1. Introduction
Population models are important tools in life sciences. Several types of them are commonly employed in the literature, and population balance models are one of the most applied models. These models describe the evolution of the population by considering that it is structured by means of several physiological characteristics. Their main features can be found in [1–3] and references therein. Cell populations can also be described in this framework. They are structured, inter alia, by the age spent in the cell cycle, cell size, or other features such as the content of groups of proteins called cyclin and intensity of certain markers. For a reference, we can mention recent literature which describes very different problems such as oscillations in a cyclin content structured model [4], the growth of yeast populations in morphologically structured ones [5, 6], or gene expression in a label-structured population [7]. References therein also provide information about the study of cell populations by means of structured models. We consider the linear size-structured population model proposed by Diekmann et al. [8]. This is a starting point in the study of more complex problems. In this model, cell mitosis happens in a symmetric way and a cell does not divide until it reaches a minimal size . This means that there is a positive minimum cell size. However, we consider that there must be a maximal size, normalized to , at which point every cell might divide or die. The model we study is given by a conservation law , , a boundary condition and an initial size distribution The independent variables and represent size and time, respectively. The dependent variable is the size-specific density of cells with size at time and we assume that the size of any individual varies according to the following ordinary differential equation: The nonnegative functions , , and represent the growth, mortality, and division rate, respectively. These are usually called the vital functions and define the life history of the individuals. Note that all of them depend on the size (the internal structuring variable). We should point out that, in (1), the reproduction process into two equal parts has been considered in the two terms in which the division rate appears. Here, we should note that the term is interpreted as zero whenever . We perform this feature with the use of functions and extended with the value zero on the interval . Condition (2) reflects that cells with a size less than cannot exist and is a consequence of the fact that cells only divide after the minimal size .
In accordance with accepted biological wisdom, there exists a maximum size. This means that cells will divide or die with probability one before reaching it. Thus, if we consider the survival property, that is, the probability that an individual of size reaches size , the hypothesis of considering a maximum size implies that . One of the forms in which this fact could be reflected consists of taking into account the growth functions introduced by Von Bertalanffy. These kinds of functions satisfy , which is enough to verify the required condition whenever and are positive and bounded functions. Note that if is a continuous function defined in , then this hypothesis implies that . Moreover, the solution to the problem must satisfy , , because we suppose that initially there are no cells of maximum size [2].
In general, physiologically structured population models are difficult to solve. Although theoretical properties of the models such as existence, uniqueness, smoothness of solutions, and long-time behaviour (with the study of steady states and their stability) could be studied without a solution expression, the knowledge of their qualitative or quantitative behaviour in a more tangible way is sometimes necessary. Therefore, numerical methods provide a valuable tool to obtain such information. In the case of general structured population models, many numerical methods have been proposed to solve them (see [9, 10] and references therein), and the difficulties found in their convergence analysis can be observed in [11], for instance.
In the case of population balance models such as (1)–(3), Liou et al. [12] proposed an alternative procedure for their solution based on a successive generations approach that provides analytical solutions in some cases. However, numerical integration may be necessary for more complicated situations. Mantzaris et al. [13] presented a finite difference scheme and Angulo and López-Marcos [14] a characteristic curve scheme with the first convergence analysis. However, until that moment, a maximum size for the cells was not considered. For the case of a bounded size interval, the works of Mantzaris et al. [15–17] provided a broad comparison of numerical methods based on finite differences and spectral and finite elements method, respectively. In their work, the authors compared the efficiency of the methods presented but they did not demonstrate their convergence and did not pay attention to either the compatibility of the initial and boundary conditions or the discontinuities caused by the maximum size. Note that the lack of smoothness properties of the solution would negatively affect the efficiency of such higher order methods. In our previous work [18], we formulated two first-order numerical procedures, a finite difference scheme and a characteristics method, and analyzed completely their convergence. The work supplied a detailed tracing of the different discontinuities arising in the simulation.
When selecting a numerical method, efficiency must be taken into account [19]. In general, on a long-time integration (e.g., see a study of the stable size distribution in [18]), the use of methods that preserve some of the qualitative properties of the solution can perform better and, thus, characteristic curves methods would be good candidates. In this paper, we consider a novel characteristics method based on the discretization of the integral representation of the solution to the problem along the characteristics lines. This procedure was previously used in [18] for this problem, obtaining a valuable first-order method. Here, in order to produce a second-order scheme, we consider a different discretization of the integral representation to the solution. Second-order methods maintain good compromise between the required smoothness of the vital functions based on realistic biological data and the efficiency of the numerical schemes. Nevertheless, this alternative discretization produces an implicit numerical method that, in principle, increases the stability property but also the computational cost. However, due to the special structure of the problem, a suitable implementation of the numerical method provides a cheaper explicit procedure.
In Section 2, we describe the numerical method to approximate the solution to (1)–(3) and comment upon the efficient implementation being used. Section 3 is devoted to the convergence analysis of the method. In Section 4, we present some numerical experiments which confirm the theoretical results and describe the performance of the numerical method in different situations related to the lack of smoothness properties of the vital functions and the solution.
2. Numerical Method
As we mentioned in Section 1, there are some schemes proposed to obtain the solution to (1)–(3). Most of them are of first-order convergence. On the one hand, this convergence property produces a lack of efficiency which can be reduced with higher order methods. On the other hand, the smoothness of the solution to (1)–(3) is not as high as these last schemes demand. However, second-order methods present a good balance: they enhance the efficiency even with a lack of regular data.
Here, we introduce an overall second-order numerical method which integrates the problem along the characteristic curves. It employs a theoretical representation of the solution to (1)–(3) whose framework was introduced in [18]. Therefore, following such work, we define and denote by the characteristic curve which takes the value at the time instant of (1). It is the solution to the following initial value problem: In this way, the solution to (1) is given by Note that, in this new layout, we have to solve two types of problems: the integration of the equation that defines the characteristic curves (6) and the solution to (7) which provides the solution to the problem along the characteristics. We use discretization procedures in order to solve them.
We consider the numerical integration of model (1)–(3) along the time interval . Thus, given a positive integer , we define and introduce the discrete time levels , . We begin with the integration of (6) which provides the grid on the space variable (size) of the method. This grid is nonuniform and invariant with time, because the growth rate function is, explicitly, independent of the time variable. However, note that it depends on time implicitly conditioned on cell size. It is usually called the natural grid [9]. In this work, we approximate such a grid by using a second-order scheme for the numerical integration of (6). More precisely, the modified Euler method provides the following approximation to the natural grid: Integer represents the index of the last grid point computed at the size interval and is chosen to satisfy the condition , with and being suitable constants (we refer to [9] for further details). Note that the points and , , , belong to the same numerical characteristic curve. Finally, we fix the last grid point .
Then, denoting , , , let be a numerical approximation to . We propose a one-step method in order to obtain it. Therefore, starting from an approximation to the initial data (3) of the problem, for example, the grid restriction of the function , the numerical solution at a new time level is described in terms of the previous one. Such a general step is obtained by means of the following second-order discretization of (7): the integrals are approached by the trapezoidal quadrature rule. For , In the previous expression, and represent approximations to the solutions at sizes and (not included in the discrete grid) and times and , respectively. So, in order to keep the second order, we compute them by linear interpolation based on the nearest grid points. More precisely, for the computation of , approximation to the solution at , and time , first we look for the index so that . Thus, Obviously, the approximating values at the minimum and maximum sizes are The numerical procedure seems to be implicit. However, if we compute the approximations at the new time level downwards (i.e., first using (11), then from to using (9), and finally using (11)), it results in an explicit procedure. The reason is that the right hand side values in (9) corresponding to the time are either zero or previously computed.
3. Convergence Analysis
In this section, we carry out the convergence analysis of the scheme. It is based on the properties of consistency and stability of the method. Henceforth, will denote a positive constant which is independent of , () and (); possibly has different values in different places.
The local discretization error is given by the following equation:
Lemma 1. Let be three times continuously differentiable; let functions , , and be two times continuously differentiable. Thus, as , the following estimates hold:
Proof. From (12) and (7), by adding and subtracting suitable terms in the expression of the local discretization error, we obtain the following bound: , . Now, taking into account the convergence properties of the trapezoidal quadrature rule, the modified Euler method, the interpolation procedure, and the assumed regularity of functions , , , and , we obtain the estimate (remember that and must be considered smooth in their extended domains).
In the following result, we prove the convergence of the numerical method. We denote the error produced by the numerical approximation by , where are the nodal values of the theoretical solution and are the numerical approximations obtained by means of the numerical method. And , .
Theorem 2. Under the hypotheses of Lemma 1, if , as , then
Proof. From (9) and (12), we have , . Taking into account the smoothness properties of the functions and , we arrive at , , and , . Then, by means of a recursive argument, we obtain , . Therefore, . By the discrete Gronwall lemma and using (13), we arrive at , and the estimate holds.
4. Numerical Results
Now, we present some different experiments in order to test the numerical method. We consider the minimum size at which a cell divides as . We suppose that there is no cellular death. Therefore, , and we choose the size-specific growth rate as .
Test Problem 1. In the first experiment, we take the size-specific division rate function We consider that the function vanishes out of interval . Coefficient is chosen in order to ensure that the maximum value of is . Note that the extended function is two times continuously differentiable.
In order to avoid discontinuities caused by an incompatible initial condition, we take satisfying . In this first experiment, we opt for with the corresponding extension out of the interval . Coefficient is chosen in order to ensure that the maximum value of is .
We do not know the analytical solution to the problem, so, in order to compare, we take as an exact solution the computed approximation with a sufficiently small value of the size step . In the experiment, we compute the solution at the final time with . If we analyze the (approximated) second derivative of the computed solution, we observe the required regularity in the hypotheses of the convergence result (see Figure 1).
In Table 1, we present the results obtained with the method for different values of the step size. For each , we compare at the final time the numerical solution being computed, , with the representation of the numerical solution corresponding to at the coarsest grid obtained with , . The second column in Table 1 shows the maximum error at the different discrete sizes; that is, The third column shows the numerical order of convergence, which we compute with the formula
Results in Table 1 clearly confirm the expected second-order of convergence.
Test Problem 2. In the second experiment, we take a size-specific division rate function not satisfying the smoothness required in the assumptions of the convergence result. For example, Again, we consider that vanishes out of the interval . Coefficient is chosen in order to ensure that the maximum value of is . Note that, now, this extended function is discontinuous.
Assuming the same initial data (25) as in the previous experiment, we obtain the results presented in Table 2. Again, we observe a second-order convergence.
Perhaps these numerical results can be explained by taking into account that the (approximated) second derivative of the product is continuous as can be observed in Figure 2. A detailed revision of the convergence result (which we do not include for the sake of simplicity) allows us to ensure the same estimative by relaxing the hypothesis, assuming that the extension of the product is two times continuously differentiable instead of the extension of .
Test Problem 3. We take the size-specific division rate function (28) of the previous experiment not satisfying the smoothness assumptions of the convergence result. But now, we consider the following initial data (producing a nonsufficiently smooth solution, as we will show), and the corresponding extended function. Coefficient is chosen in order to ensure that the maximum value of is . Note that, now, the extended function is not differentiable.
In this case, at , we observe that the (approximated) second derivative of is discontinuous (Figure 3(a)). And we see the same behavior for the (approximated) second derivative of the product (Figure 3(b)). Now, the previous convergence analysis is not valid.
(a)
(b)
Table 3 shows the results obtained with the method. We conclude that there is still convergence in the numerical approximation. However, we do not observe a well-defined order.
5. Conclusions
The study of cell populations by means of the use of size-structured models is a current topic. Its numerical integration exhibits a great development in obtaining qualitative or quantitative information about the solution. However, there is a lack of attention to certain important problems which take part in such integration. First, it is necessary to respect the individual maximum size, which is biological wisdom. Second, it is appropriate to increase the efficiency of the integration with the use of higher order methods. However, it is difficult to blend these two components due to the lack of smoothness which could appear in general situations.
We have proposed a new numerical method to attain the solution to (1)–(3). We have proved its second-order convergence, and we have corroborated it experimentally. We also observe, numerically, its robustness in different situations.
The numerical method proposed in this work can be extended to different situations such as the uneven division case [20], assuming density independent division functions. In such circumstances, the birth term in the equation involves an integral that must be approximated by means of a suitable quadrature rule. The resulting method differs from the scheme considered in this paper. For example, the discretization of the nonlocal term in the equal fission case requires an interpolation procedure, but it is not necessary in the uneven case. The numerical scheme designed comes from an implicit discretization but, again, a downward implementation provides an explicit numerical procedure. However, for this more general model, the convergence proof of the numerical method could be harder than that for the equal partitioning case.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the reviewers for their constructive and helpful suggestions. This work was supported in part by Projects MTM2011-25238 and MTM2014-56022-C2-2-P of the Ministerio de Economía y Competitividad (Spain) and VA191U13 of the Junta de Castilla y León (Spain).