Optimal control model of tumor treatment in the context of cancer stem cell

1.   Introduction

Coagulation-fragmentation (CF) equations have been used to model many physical and biological phenomena ^[1,2]. In particular, when combined with transport terms, these equations can be used to model the population dynamics of oceanic phytoplankton ^[3,4,5]. Setting such models in the space of Radon measures allows for the unified study of both discrete and continuous structures. Not only are the classical discrete and continuous CF equations special cases of the measure valued model (as shown in ^[6]), but this setting allows for a mixing of the two structures, which has become of interest in particular applications ^[7,8].

With the above applications in mind, numerical schemes to solve CF equations are of great importance to researchers. In particular, finite difference methods offer numerical schemes which are easy to implement and approximate the solution with a high order of accuracy. The latter benefit is especially important in the study of stability and optimal control of such equations.

The purpose of this article is to make improvements on two of the three first-order schemes presented in ^[9], namely the fully explicit and semi-implicit schemes. These schemes are shown to have certain advantages and disadvantages as discussed in the aforementioned study. In particular, the fully explicit scheme has the qualitative property of conservation of mass through coagulation. On the other hand, the semi-implicit scheme has a more relaxed Courant–Friedrichs–Lewy (CFL) condition, which does not depend on the initial condition. We have decided not to attempt to improve the third scheme presented in ^[9] as there does not seem to be a significant advantage of the named conservation law scheme to outweigh the drastic computational cost. The main improvement here is to lift-up these two first-order schemes to second-order ones on the space of Radon measures; however, as this state space contains singular elements (including point measures), the improvement of these schemes must be handled with care. As shown in ^[10], discontinuities and singularities in the solution can cause drastic changes in not only the order of convergence of the scheme, but also in the behavior of the scheme. To address these issues, we turn to a high resolution scheme studied with classical structured population models (i.e. without coagulation-fragmentation) in ^[11,12,13]. This scheme makes use of a minmod flux limiter to control any oscillatory behavior of the scheme caused by irregularities. With this new flux, we show that it is possible for second order convergence rates to be obtained for continuous density solutions. However, as the solutions become more irregular, one should expect the convergence rate to decline. Such a phenomenon is demonstrated in ^[10,13], and we direct the reader to these manuscripts for more discussion.

The layout of the paper is as follows. In Section 2, we present the notation and preliminary results about the model and state space used throughout the paper. In Section 3, we describe the model and state all assumptions imposed on the model parameters. In Section 4, we present the numerical schemes, their CFL conditions, and state the main theorem of the paper. In Section 5, we test the convergence rate of the schemes against well-known examples. In Section 6, we provide a conclusion and in the Appendix (Section 7) we provide proofs for some of our results.

2.   Notation

We make use of standard notations for function spaces. The most common examples of these are $C^1(\mathbb{R}^+)$ for the space of real valued continuously differentable functions and $W^{1, \infty}(\mathbb{R}^+)$ for the usual Sobelov space. The space of Radon measures will be denoted with $\mathcal{M}(\mathbb{R}^+)$, with $\mathcal{M}^+(\mathbb{R}^+)$ representing its positive cone. This space will be equipped with the Bounded-Lipschitz (BL) norm given by

Another norm of interest to this space is the well studied Total Variation (TV) norm given by

For more information about these particular norms and their relationship we direct the reader to ^[14,15]. For lucidity, we use operator notation in place of integration when we believe it necessary, namely

where the set $A$ is the support of the measure $\mu$. Finally, we denote the minmod function by $\text{mm}(a, b)$ and use the following definition

3.   Model and assumptions

The model of interest is the size-structured coagulation fragmentation model given by

where $\mu(t) \in \mathcal{M}^+(\mathbb{R}^+)$ represents individuals' size distribution at time $t$ and the functions $g, d, \beta$ are their growth, death, and reproduction rate, respectively. The coagulation and fragmentation processes of a population distributed according to $\mu \in \mathcal{M}^+(\mathbb{R}^+)$ are modeled by the measures $K[\mu]$ and $F[\mu]$ given

and

for any test function $\phi$. Here, $\kappa(x, y)$ is the rate at which individuals of size $x$ coalesce with individuals of size $y$, $a(y)$ is the global fragmentation rate of individuals of size $y$, and $b(y, \cdot)$ is a measure supported on $[0, y]$ such that $b(y, A)$ represents the probability a particle of size $y$ fragments to a particle with size in the Borel set $A$.

Definition 3.1. Given $T\geq 0$, we say a function $\mu \in C([0, T], \mathcal{M^+}(\mathbb{R}^+))$ is a weak solution to (3.1) if for all $\phi \in (C^1 \cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$ and for all $t\in [0, T]$, the following holds:

For the numerical scheme, we will restrict ourselves to a finite domain, $[0, x_{\max}]$. Thus, we impose the following assumptions on the growth, death and birth functions:

$({\rm{A1}})$ For any $R > 0$, there exists $L_R > 0$ such that for all $\| \mu_i\|_{TV}\leq R$ and $t_i\in [0, \infty)$ ($i = 1, 2$) the following hold for $f = g, d, \beta,$

$({\rm{A2}})$ There exists $\zeta > 0$ such that for all $T > 0$,

$({\rm{A3}})$ For all $(t, \mu)\in [0, \infty)\times \mathcal{M^+}(\mathbb{R}^+)$,

for some large $x_{\max} > 0$.

We assume that the coagulation kernel $\kappa$ satisfies the following assumption:

$({\rm{K1}})$ $\kappa$ is symmetric, nonnegative, bounded by a constant $C_\kappa$, and globally Lipschitz with Lipschitz constant $L_\kappa$.

$({\rm{K2}})$ \(\kappa(x, y) = 0 \text{ whenever } x+y > x_{\max}.\)

We assume that the fragmentation kernel satisfies the following assumptions:

$({\rm{F1}})$ $a\in W^{1, \infty}(\mathbb{R}^+)$ is non-negative,

$({\rm{F2}})$ for any $y\ge 0$, $b(y, dx)$ is a measure such that

$({\rm{i}})$ $b(y, dx)$ is non-negative and supported in $[0, y]$, and there exist a $C_b > 0$ such that $b(y, \mathbb{R}^+)\le C_b$ for all $y > 0$,

$({\rm{ii}})$ there exists $L_b$ such that for any $y, \bar y\ge 0$,

$({\rm{iii}})$ for any $y\ge 0$,

The existence and uniqueness of mass conserving solutions of model (3.1) under these assumptions were established in ^[6].

4.   Numerical methods

We adopt the numerical discretization presented in ^[6]. For some fixed mesh sizes $\Delta x, \Delta t > 0$, we discretize the size domain $[0, x_{\max}]$ with the cells

and

We denote the midpoints of these grids by $x_j$. The initial condition $\mu_0 \in \mathcal{M}^+(\mathbb{R}^+)$ will be approximated by a combination of Dirac measures

We first approximate the model coefficients $\kappa$, $a$, $b$ as follows. For the physical ingredients, we define

for $i, j\ge 1$, and

(with the natural modifications for $\kappa^{\Delta x}_{0, j}$ and $\kappa^{\Delta x}_{i, 0}$, $i\ge 1$). We then let $a^{\Delta x}\in W^{1, \infty}(\mathbb{R}^+)$ and $\kappa^{\Delta x}\in W^{1, \infty}(\mathbb{R}^+\times \mathbb{R}^+)$ be the linear interpolation of the $a^{\Delta x}_i$ and $\kappa^{\Delta x}_{i, j}$, respectively. Finally, we define the measure $b^{\Delta x}(x_j, \cdot)\in \mathcal{M}^+(\Delta x \mathbb{N})$ by

and then $b^{\Delta x}(x, \cdot)\in \mathcal{M}^+(\Delta x \mathbb{N}_0)$ for $x\ge 0$ as the linear interpolate between the $b^{\Delta x}(x_j, \cdot)$. When the context is clear, we omit the $\Delta x$ from the notation above.

We make use of these approximations to combine the high-resolution scheme presented in ^[13] with the fully explicit and semi-implicit schemes presented in ^[9]. Together, these schemes give us the numerical scheme

where the flux term is given by

the fragmentation term, $\mathcal{F}_{j, k}$, is given by

and the coagulation term, $\mathcal{C}_j$, is either given by an explicit discretization as

or by an implicit one as

As discussed in ^[9], the explicit scheme which uses (4.4) to approximate the coagulation term and the semi-implicit scheme which instead uses (4.5) to approximate the coagulation term behave differently with respect to the mass conservation and have different Courant–Friedrichs–Lewy (CFL) conditions. The assumed CFL condition for the schemes are

where $\bar{\zeta} = \max\{ \zeta, \|a\|_{W^{1, \infty}}\}$, $C_a = \|a\|_\infty$. The CFL conditions above are similar to those used in ^[9], but are adjusted due to the flux limiter term as in ^[13]. It is clear that the semi-implicit scheme has a less restrictive and simpler CFL condition than the explicit scheme. In particular, the CFL condition of the semi-implicit scheme is independent on the initial condition, unlike its counterpart. The trade off for this is a loss of qualitative behavior of the scheme in the sense of mass conservation. Indeed as shown in ^[9], when $\beta = d = g = a = 0$, the semi-implicit coagulation term does not conserve mass, whereas the explicit term does. As shown in the appendix, this loss is controlled by the time step size, $\Delta t$.

It is useful to define the following coefficients:

and

Notice, $|A_j^k| \leq \frac{3\Delta t}{2\Delta x}\zeta$ and $A_j^k -B_j^k \geq 0$ as

Scheme (4.1) can then be rewritten as

Depending on the choice of coagulation term, this formulation leads to either

for the explicit term, $\mathcal{C}^{\text{exp}}_{j, k}$, or

for the implicit term, $\mathcal{C}^{\text{imp}}_{j, k}$.

For these, schemes, we have the following Lemmas which are proven in the appendix:

Lemma 4.1. For each $k = 1, 2, \dots, \bar{k}$,

● $m_j^k \geq 0$ for all $j = 1, 2, \dots J$,

● $\| \mu^k_{\Delta x}\|_{TV} \leq \|\mu_0 \|_{TV} \exp((\zeta + C_b C_a)T).$

Lemma 4.2. For any $l, p = 1, 2, \dots, \bar{k}$,

Using the above two Lemmas, we can arrive at analogous results for the linear interpolation (4.10):

Thus by the well know Ascoli-Arzela Theorem, we have the existence of a convergent subsequence of the net $\{\mu_{\Delta x}^{\Delta t}(t) \}$ in $C([0, T], \mathcal{M}^+([0, x_{\max}])$. We now need only show any convergent subsequence converges to the unique solution (3.2).

Theorem 4.1. As $\Delta x, \Delta t\to 0$ the sequence $\mu_{\Delta x}^{\Delta t}$ converges in $C([0, T], \mathcal{M^+}([0, x_{\max}]))$ to the solution of (3.1).

Proof. By multiplying (4.1) by a superfluously smooth test function $\phi \in (W^{1, \infty} \cap C^2)([0, T]\times \mathbb{R})$, denoting $\phi_j^k : = \phi (k\Delta t, x_j)$, summing over all $j$ and $k$, and rearranging we arrive at

The left-hand side of equation (4.11) was shown in ^[13] to be equivalent to

where $o(1) \longrightarrow 0$ as $\Delta t, \Delta x \longrightarrow 0$.

The right-hand side of (4.11) was shown in ^[9] to be equal to

Making use of results, it is then easy to see (4.11) is equivalent to

Passing the limit as $\Delta t, \Delta x \longrightarrow 0$ along a converging subsequence, we then obtain that equation (3.2) holds for any $\phi \in (C^2\cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$ with compact support. A standard density argument shows that equation (3.2) holds for any $\phi \in (C^1\cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$. As the weak solution is unique ^[6], we conclude the net $\{\mu_{\Delta x}^{\Delta t}\}$ converges to the solution of model (3.1). □

We point out that while these schemes are higher-order in space, they are only first order in time. To lift these schemes into a second-order in time as well, we make use of the second-order Runge-Kutta time discretization ^[16] for the explicit scheme and second-order Richardson extrapolation ^[17] for the semi-implicit scheme.

5.   Numerical examples

In this section, we provide numerical simulations which test the order of the explicit and semi-implicit schemes developed in the previous sections. We test each component separately, beginning first with a pure coagulation equation in example 1 (where we set $g = \beta = d = a = 0$), then a pure fragmentation equation in example 2 (where we set $g = \beta = d = \kappa = 0$). In example 3, we consider all components of model (3.1) including the boundary term which is implemented as in scheme (4.7). For readers interested in the schemes performance in the absence of the coagulation-fragmentation processes, we direct you to ^[11,12,13]. For each example, we give the BL error and the order of convergence. To appreciate the gain in the order of convergence compared to those studied in ^[9], which are based on a first order approximation of the transport term, we add some of the numerical results from the scheme presented in ^[9].

In some of the following examples, the exact solution of the model problem is given. In these cases, we approximate the order of accuracy, $q$, with the standard calculation:

where $\mu$ represents the exact solution of the examples considered. In the cases where the exact solutions are unknown, we approximate the order by

and we report the numerator of the log argument as the error. The metric $\rho$ we use here was introduced in ^[18] and is equivalent to the BL metric, namely

for some constant $C$ (dependent on the finite domain). As discussed in ^[18], this metric is more efficient to compute than the BL norm and maintains the same order of convergence. An alternative to this algorithm would be to make use of the algorithms presented in ^[19], where convergence in the Fortet-Mourier distance is considered.

Example 1 In this example, we test the quality of the finite difference schemes against coagulation equations. To this end, we take $\kappa (x, y) \equiv 1$ and $\mu_0 = e^{-x} dx$ with all other ingredients set to $0$. This example has an exact solution given by

see ^[20] for more details. The simulation is performed over the truncated domain $(x, t) \in [0, 20]\times [0, 0.5]$. We present the BL error and the numerical order of convergence for both schemes in Table 1.

Example 2 In this example, we test the quality of the finite difference scheme against fragmentation equations. We point out that in this case, the two schemes are identical in the spacial component. For this demonstration, we take $\mu_0 = e^{-x}dx$, $b(y, \cdot) = \frac{2}{y}dx \text{ and } a(x) = x$. As given in ^[21], this problem has an exact solution of

The simulation is performed over the finite domain $(x, t) \in [0, 20]\times [0, 0.5]$. We present the BL error and the numerical order of convergence for both schemes in Table 2. Note as compared to coagulation, the fragmentation process is more affected by the truncation of the domain. This results in the numerical order of the scheme being further from 2 than example 1.

Example 3 In this example, we test the schemes against the complete model (i.e., with all biological and physical processes). To this end, we take $\mu_0 = e^{-x}dx$, $g(x) = 2-2e^{x-20}$, $\beta (x) = 2$, $d(x) = 1$, $\kappa (x, y) = 1$, $a(x) = x$, and $b(y, \cdot) = \frac{2}{y}$. The simulation is performed over the finite domain $(x, t) \in [0, 20]\times [0, 0.5]$. To our knowledge, the solution of this problem is unknown.

Example 4 As mentioned in ^[9], the mixed discrete and continuous fragmentation model studied in ^[7,8], with adjusted assumptions, is a special case of model (3.1). Indeed, by removing the biological and coagulation terms and letting the kernel

with $\text{supp} \ b^c(y, \cdot) \subset [Nh, y]$ for some $h > 0$, we have the mixed model in question. We wish to demonstrate the finite difference scheme presented here maintains this mixed structure.

To this end, we take the fragmentation kernel

with initial condition $\mu = \sum_{i = 1}^5 \delta_i + \chi_{[5,15]}(x)dx$, where $\chi_A$ represents the characteristic function over the set $A$. This is similar to some examples in ^[8], where more detail and analysis are provided. In Figure 1, we present the simulation of this example. Notice, the mixed structure is preserved in finite time. For examples of this type, the scheme could be improved upon by the inclusion of mass conservative fragmentation terms similar to those presented in ^[6].

6.   Conclusion

In this paper, we have lifted two of the first order finite difference schemes presented in ^[9] to second order high resolution schemes using flux limiter methods. The difference between both schemes is only found in the coagulation term, where the semi-implicit scheme is made linear. In context of standard structured population models (i.e. without coagulation or fragmentation), these type of schemes have been shown to be well-behaved in the presences of discontinuities and singularities. This quality makes them a well fit tool for studying PDEs in spaces of measures. We prove the convergence of both schemes under the assumption of natural CFL conditions. The order of convergence of both schemes is then tested numerically with previously used examples.

In summary, the schemes preform as expected in the presence of smooth initial conditions. In all such simulations, the numerical schemes presented demonstrate a convergence rate of order 2. For simulations with biological terms, this convergence rate is expected to drop when singularities and discontinuities occur, as demonstrated in ^[13]. Mass conservation of the schemes, an important property for coagulation/fragmentation processes, is discussed in detail in ^[6,9].

Acknowledgments

The research of ASA is supported in part by funds from R.P. Authement Eminent Scholar and Endowed Chair in Computational Mathematics at the University of Louisiana at Lafayette. RL is grateful for the support of the Carl Tryggers Stiftelse via the grant CTS 21:1656.

7.   Appendix

7.1.   Proof of Lemmas 4.1 and 4.2

In this section, we present the proofs of Lemmas 4.1 and 4.2 for the explicit coagulation term. The semi-implicit term follows from similar arguments in the same fashion as ^[9].

Proof of Lemma 4.1

Proof. We first prove via induction that for any $k = 1, 2, \dots, \bar{k},$ $\mu_{\Delta x}^k$ satisfies the following:

$({\rm{i}})$ $\mu^k_{\Delta x} \in \mathcal{M}^+(\mathbb{R}^+)$ i.e. $m_j^k \geq 0$ for all $j = 1, \dots, J$,

$({\rm{ii}})$ $\|\mu_{\Delta x}^k\|_{TV} \leq \|\mu^0_{\Delta x}\|_{TV} (1+(\zeta + C_b C_a)\Delta t)^k.$

Then, the TV bound in the Lemma follows from standard arguments (see e.g. Lemma 4.1 in ^[9]). We prove this Theorem for the choice of the explicit coagulation term, $\mathcal{C}_{j, k}^{\text{exp}},$ as the implicit case is similar and more straight forward.

We begin by showing that $m_j^{k+1} \geq 0$ for every $j = 1, 2, \dots, J$. Notice by way of (4.8), this reduces down to showing

Indeed, by the CFL condition (4.6), induction hypothesis, and

we arrive at the result.

For the TV bound, we have since the $m_j^k$ are non-negative, $\|\mu_{\Delta x}^k\|_{TV} = \sum_{j = 1}^{J} m_j^k$. By rearranging (4.8) and summing over $j = 1, 2, \dots, J$ we have

To bound the right-hand side of equation (7.1), we directly follow the arguments of Lemma 4.1 in ^[9] which yields

Using the induction hypothesis, we obtain $\|\mu^{k+1}_{\Delta x}\|_{TV}\le \|\mu^0_{\Delta x}\|_{TV} (1+(\zeta + C_b C_a)\Delta t)^{k+1}$ as desired. □

Proof of Lemma 4.2

Proof. For $\phi \in W^{1, \infty}(\mathbb{R}^+)$ with $\|\phi\|_{W^{1, \infty}} \leq 1$, and denoting $\phi_j: = \phi(x_j)$, we have for any $k$,

Let $C$ be the right-hand side of the TV-bound from Lemma 4.1, we then see

Moreover, since $g_J^k = 0$ the sum in the right-hand side takes the form

We thus obtain

Taking the supremum over $\phi$ gives $\|\mu^{k+1}_{\Delta x}-\mu^k_{\Delta x}\|_{BL}\le L\Delta t$ for any $k$. The result follows. □

7.2.   Estimate of mass loss for the implicit coagulation discretization given by ${\mathcal C}^{\text{imp}}_{j, k}$

In this section, we consider the semi-implicit scheme (4.9) without any biological ingredients or fragmentation (i.e., $a, g, d, \beta = 0$) where mass is not conserved (observed by the numerical experiments in Section 5) and show this change is controlled by the time step. For a bound on the loss of mass via fragmentation, we direct the reader to Section 6.1 of ^[6]. Multiplying (4.9) by $x_j$ and summing over $j$, we can arrive at

In [6, Section 6.1] it was shown that the explicit scheme conserves mass through the coagulation process, i.e.,

Adding this to the previous equation we have

where in the last equality, we change the order of integration and introduce the new index $l = j-i$. Noticing that due to the uniform mesh size $x_{l+i} = x_l + x_i$, we can split the second term on the right-hand side and obtain the equation

Since $x_j \leq x_{\max}$, we can bound the last term on the right-hand side

Using Lemmas 4.1 and 4.2, we have the estimate