# A Second-Order Method for the Numerical Integration of a Size-Structured Cell Population Model

Authors: O. Angulo; J. C. López-Marcos; M. A. López-Marcos
Journal: Abstract and Applied Analysis
 (2015)

---

## 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.


---

## Body

## 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 a>0. This means that there is a positive minimum cell size. However, we consider that there must be a maximal size, normalized to x=1, at which point every cell might divide or die. The model we study is given by a conservation law (1)utx,t+gxux,tx=-μ(x)u(x,t)-b(x)u(x,t)+4b(2x)u(2x,t),a/2<x<1, t>0, a boundary condition (2)ua2,t=0,t>0,and an initial size distribution (3)u(x,0)=φ(x),a2≤x≤1.The independent variables x and t represent size and time, respectively. The dependent variable u(x,t) is the size-specific density of cells with size x at time t and we assume that the size of any individual varies according to the following ordinary differential equation: (4)dxdt=g(x).The nonnegative functions g, μ, and b 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 x (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 4b(2x)u(2x,t) is interpreted as zero whenever x>1/2. We perform this feature with the use of functions u and b extended with the value zero on the interval [1,2]. Condition (2) reflects that cells with a size less than a/2 cannot exist and is a consequence of the fact that cells only divide after the minimal size a>0.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 sizex0 reaches size x, (5)Π(x)=exp⁡-∫x0xμ(s)+b(s)g(s)ds,a2≤x0≤x<1,the hypothesis of considering a maximum size implies that limx→1-Π(x)=0. 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 ∫x01(ds/gs)=∞, which is enough to verify the required condition whenever μ and b are positive and bounded functions. Note that if g is a continuous function defined in a/2,1, then this hypothesis implies that g(1)=0. Moreover, the solution to the problem must satisfy u(1,t)=0, t>0, 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 Section2, 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 Section1, 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 μ∗(x)=g′(x)+μ(x)+b(x) and denote by x(t;t∗,x∗) the characteristic curve which takes the value x∗ at the time instant t∗ of (1). It is the solution to the following initial value problem: (6)ddtx(t;t∗,x∗)=g(x(t;t∗,x∗)),t>t∗,x(t∗;t∗,x∗)=x∗.In this way, the solution to (1) is given by (7)uxt;t∗,x∗,t=u(x∗,t∗)·exp⁡-∫t∗tμ∗xτ;t∗,x∗dτ+∫t∗texp⁡-∫τtμ∗xs;t∗,x∗dsmmm·4b(2xτ;t∗,x∗)mmm·u(2xτ;t∗,x∗,τ)dτ,t≥t∗.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 [0,T]. Thus, given a positive integer N, we define k=T/N and introduce the discrete time levels tn=nk, 0≤n≤N. 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 thenatural 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: (8)x0=a2,xj+1=xj+k2gxj+gxj+kg(xj),hhhhhhhhhhhhhhhhhhh0≤j≤J-1.Integer J represents the index of the last grid point computed at the size interval and is chosen to satisfy the condition K0k≤1-xJ≤K1k, with K0 and K1 being suitable constants (we refer to [9] for further details). Note that the points (xj,tn) and (xj+1,tn+1), 0≤j≤J-1, 0≤n≤N-1, belong to the same numerical characteristic curve. Finally, we fix the last grid point xJ+1=1.Then, denotingujn=u(xj,tn), 0≤j≤J+1, 0≤n≤N, let Ujn be a numerical approximation to ujn. 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 0≤n≤N-1, (9)Uj+1n+1=exp⁡-k2μ∗xj+μ∗xj+1·Ujn+2kb(2xj)U-2·jn+2kb(2xj+1)U-2·(j+1)n+1,j=0,…,J-1.In the previous expression, U-2·jn and U-2·(j+1)n+1 represent approximations to the solutions at sizes 2xj and 2xj+1 (not included in the discrete grid) and times tn and tn+1, 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 U-2·lm, approximation to the solution at 2xl, and time tm, first we look for the index M so that xM-1<2xl≤xM. Thus, (10)U-2·lm=UM-1m+UMm-UM-1mxM-xM-12xl-xM-1,if2xl<1,0,if2xl≥1.Obviously, the approximating values at the minimum and maximum sizes are (11)U0n+1=UJ+1n+1=0.The numerical procedure seems to be implicit. However, if we compute the approximations at the new time level tn+1 downwards (i.e., first UJ+1n+1 using (11), then Uj+1n+1 from J-1 to 0 using (9), and finally U0n+1 using (11)), it results in an explicit procedure. The reason is that the right hand side values in (9) corresponding to the time tn+1 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,C will denote a positive constant which is independent of k, n (0≤n≤N) and j (0≤j≤J+1); C possibly has different values in different places.The local discretization error is given by the following equation: (12)τj+1n+1=1kuj+1n+1-exp⁡-k2μ∗xj+μ∗xj+1mmmm-k2μ∗xj+μ∗xj+1·ujn+2kb(2xj)u-2·jn-2kb(2xj+1)u-2·(j+1)n+1,mmmmmmmmmmmmc0≤j≤J-1,0≤n≤N-1.Lemma 1.
Letg be three times continuously differentiable; let functions μ, b, and u be two times continuously differentiable. Thus, as k→0, the following estimates hold: (13)τj+1n+1=O(k2),0≤j≤J-1,0≤n≤N-1.Proof.
From (12) and (7), by adding and subtracting suitable terms in the expression of the local discretization error, we obtain the following bound: (14)τj+1n+1≤1kuj+1n+1-uxtn+1;xj,tn,tn+1+1kujnexp⁡-∫tntn+1μ∗xτ;tn,xjdτk2mmmmmmc-exp⁡-k2μ∗xjxtn+1;tn,xjmmmmmmmmmmmmm∫tntn+1μ∗xτ;tn,xjdτk2+μ∗xtn+1;tn,xj+1kujnexp⁡-k2μ∗xj·exp⁡-k2μ∗xtn+1;tn,xjmmm-exp⁡-k2μ∗xj+1+4k∫tntn+1exp⁡-∫τtn+1μ∗xs;tn,xjdsmmmmmmm·b(2xτ;tn,xj)mmmmmmm·u(2xτ;tn,xj,τ)dτmmmcmm-k2exp⁡-∫tntn+1μ∗xs;tn,xjdsmmmcmmmmmm·b(2xj)u(2xj,tn)mmmmmcmmmm+b(2xtn+1;tn,xj)mmmmmmmmm-∫tntn+1μ∗xs;tn,xjds·u(2xtn+1;tn,xj,tn+1)+2exp⁡-∫tntn+1μ∗xs;tn,xjdsmmmm-∫tntn+1μ∗xs;tn,xjds-exp⁡-k2μ∗xj+μ∗xtn+1;tn,xj·b(2xj)u(2xj,tn)+2exp⁡-k2μ∗xj·exp⁡-k2μ∗xtn+1;tn,xjmmm-exp⁡-k2μ∗xj+1·b(2xj)u(2xj,tn)+2exp⁡-k2μ∗xj+μ∗xj+1b2xj·u(2xj,tn)-u-2·jn+b2xtn+1;tn,xj-b2xj+1·u(2xtn+1;tn,xj,tn+1)+b2xj+1·u(2xtn+1;tn,xj,tn+1)-u(2xj+1,tn+1)+b(2xj+1)u(2xj+1,tn+1)-u-2·(j+1)n+1,0≤j≤J-1, 0≤n≤N-1. 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 g, μ, b, and u, we obtain the estimate (remember that u and b 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(15)en=(e0n,…,eJn,eJ+1n),ejn=ujn-Ujn,0≤j≤J+1,0≤n≤N, where ujn are the nodal values of the theoretical solution and Ujn are the numerical approximations obtained by means of the numerical method. And (16)e-2·lm=u-2·lm-U-2·lm,xM-1<2xl≤xM≤1,0,2xl≥1,0≤l≤J, 0≤m≤N.Theorem 2.
Under the hypotheses of Lemma1, if e0∞=O(k2), as k→0, then (17)en∞≤Ck2,0≤n≤N.Proof.
From (9) and (12), we have (18)ej+1n+1=exp⁡-k2μ∗xj+μ∗xj+1ejn+2kb(2xj)e-2·jn-2kb(2xj+1)e-2·(j+1)n+1+kτj+1n+1,0≤j≤J-1, 0≤n≤N-1. Taking into account the smoothness properties of the functions μ∗ and b, we arrive at (19)ej+1n+1≤1+Ckejn+Cken∞+en+1∞+kτj+1n+1,0≤j≤J-1, 0≤n≤N-1, and (20)|ej+1n+1|≤1+Cken∞+Cken+1∞+kτn+1∞,0≤j≤J-1, 0≤n≤N-1. Then, by means of a recursive argument, we obtain (21)ejn≤1+Ckne0∞+Ck∑l=1n1+Ckn-lel∞+k∑l=1n1+Ckn-lτl∞,1≤j≤J, 1≤n≤N. Therefore, (22)en∞≤Ce0∞+k∑l=1nel∞+k∑l=1nτl∞,1≤n≤N. By the discrete Gronwall lemma and using (13), we arrive at (23)en∞≤Ce0∞+O(k2),0≤n≤N, 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 asa=1/4. We suppose that there is no cellular death. Therefore, μ(x)=0, and we choose the size-specific growth rate as g(x)=0.11-x.Test Problem 1. In the first experiment, we take the size-specific division rate function (24)b(x)=b0x-1431-x3,14≤x≤1.We consider that the function vanishes out of interval 1/4,1. Coefficient b0 is chosen in order to ensure that the maximum value of b(x) is 1. Note that the extended function is two times continuously differentiable.In order to avoid discontinuities caused by an incompatible initial condition, we takeφ satisfying φ1/8=φ′1/8=φ′′(1/8)=0. In this first experiment, we opt for (25)φ(x)=φ0x-1831-x3,18≤x≤1,with the corresponding extension out of the interval 1/8,1. Coefficient φ0 is chosen in order to ensure that the maximum value of φ(x) is 1.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 stepk. In the experiment, we compute the solution at the final time T=1 with k∗=4.8828125e-4. 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).Figure 1
Test Problem 1. Approximated second derivative ofu.In Table1, we present the results obtained with the method for different values of the step size. For each k, we compare at the final time T the numerical solution being computed, UkN, with the representation of the numerical solution corresponding to k∗ at the coarsest grid obtained with k, Uk∗N. The second column in Table 1 shows the maximum error at the different discrete sizes; that is, (26)ek=UkN-Uk∗N∞.The third column shows the numerical order of convergence, which we compute with the formula (27)s=log⁡(e2k/ek)log⁡(2).Table 1
Test Problem 1. Error and numerical convergence order.T=1.
k
Error
Order
2.5
e
-
1
2.181199
e
-
2
1.25
e
-
1
5.968912
e
-
3
1.9
6.25
e
-
2
1.484527
e
-
3
2.0
3.125
e
-
2
3.700337
e
-
4
2.0
1.5625
e
-
2
9.218746
e
-
5
2.0
7.8125
e
-
3
2.293369
e
-
5
2.0Results in Table1 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, (28)b(x)=b1x-143,14≤x≤1.Again, we consider that b(x) vanishes out of the interval 1/4,1. Coefficient b1 is chosen in order to ensure that the maximum value of b(x) is 1. 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.Table 2
Test Problem 2. Error and numerical convergence order.T=1.
k
Error
Order
2.5
e
-
1
3.975049
e
-
3
1.25
e
-
1
1.001526
e
-
3
2.0
6.25
e
-
2
2.507402
e
-
4
2.0
3.125
e
-
2
6.256562
e
-
5
2.0
1.5625
e
-
2
1.565355
e
-
5
2.0
7.8125
e
-
3
3.892902
e
-
6
2.0Perhaps these numerical results can be explained by taking into account that the (approximated) second derivative of the productbu 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 bu is two times continuously differentiable instead of the extension of b.Figure 2
Test Problem 2. Approximated second derivative ofbu.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), (29)φ(x)=φ1x-183(1-x),18≤x≤1,and the corresponding extended function. Coefficient φ1 is chosen in order to ensure that the maximum value of φ(x) is 1. Note that, now, the extended function is not differentiable.In this case, atT=1, we observe that the (approximated) second derivative of u is discontinuous (Figure 3(a)). And we see the same behavior for the (approximated) second derivative of the product bu (Figure 3(b)). Now, the previous convergence analysis is not valid.Figure 3
Test Problem 3. Approximated second derivative ofu (a) and approximated second derivative of bu (b).
(a)
(b)Table3 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.Table 3
Test Problem 3. Error and numerical convergence order.T=1.
k
Error
Order
2.5
e
-
1
3.015407
e
-
2
1.25
e
-
1
6.003792
e
-
3
2.3
6.25
e
-
2
2.142219
e
-
3
1.5
3.125
e
-
2
1.560384
e
-
4
3.8
1.5625
e
-
2
7.186373
e
-
5
1.1
7.8125
e
-
3
2.946163
e
-
5
1.3


## 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.



---

*Source: 549168-2015-08-06.xml*