A mathematical model of tumour growth with Beddington–DeAngelis functional response: a case of cancer without disease

ABSTRACT A previously published mathematical model, governing tumour growth with mixed immunotherapy and chemotherapy treatments, is modified and studied. The search time, which is assumed to be neglectable in the previously published model, is incorporated into the functional response for tumour-cell lysis by effector cells. The model exhibits bistability where a tumour-cell population threshold exists. A tumour with an initial cell population below the threshold can be controlled by the immune system and remains microscopic and asymptomatic called cancer without disease while that above the threshold grows to lethal size. Bifurcation analysis shows that (a) the chemotherapy-induced damage may cause a microscopic tumour, which would never grow to become lethal if untreated, to grow to lethal size, (b) applying chemotherapy alone requires a large dosage to be successful, (c) immunotherapy is ineffective, and (d) a combination of chemotherapy and immunotherapy can produce a synergistic effect on the outcome of a treatment.


Introduction
Cancer has become a leading cause of death worldwide. Cancer cells grow rapidly and have the ability to metastasize to other organs. The goal of cancer treatments is to destroy cancer cells without damage to healthy cells. Chemotherapy and immunotherapy are among the most common options for cancer treatment. Chemotherapy is the use of drugs to destroy cancer cells, and immunotherapy is a treatment that enhances the immune system in fighting cancer. Unfortunately, chemotherapy destroys not only the cancer cells but also the surrounding healthy cells; additionally, more effective drugs are likely to be more toxic. Immunotherapy has been shown to be ineffective in patients with solid tumours or advanced cancer [9,22,38]. Some research studies have suggested that an appropriate combination of chemotherapy and immunotherapy can produce a synergistic effect [27,29,34,38].
Mathematical modelling of tumour growth is a powerful tool to understand, predict, and improve the outcome of a treatment regimen. Many mathematical models of tumour growth use nonlinear ordinary differential equations (ODEs) [10][11][12][16][17][18][24][25][26]31,35,36]. This work is based on a mathematical model developed by de Pillis et al. [18]. The mathematical model is a system of ODEs governing the tumour growth on a cell population level with a ratio-dependent like interaction between tumour cells and cytotoxic T cells. The system also considers a combination of immunotherapy and chemotherapy treatments, and the parameter values have been estimated using experimental data from clinical research studies.
Due to the intrinsic property of the ratio-dependent functional response, the mathematical model [18] is not differentiable at the tumour-free equilibrium. Linear stability analysis fails in identifying the stability of the tumour-free equilibrium. Wei [42] has proven that the tumour-free equilibrium is stable under a range of parameter values without the treatment terms in the mathematical model. Numerical simulation has been conducted by de Pillis et al. [18] and Wei [42] for the study of cancer treatments. de Pillis et al. have performed integration of the system using different treatment regimens and concluded that a combination of chemotherapy and immunotherapy improves the efficacy and safety of the treatment because it allows the use of a small dosage of each cancer therapy. Wei [42] has performed one-dimensional bifurcation analysis and found that the unstable fixed point curve provides a good approximation of the maximum tumour burden that the immune system is able to control. When chemotherapy is used alone, the system exhibits chemotherapy-induced lymphocyte damage, which could change the stability of the tumour-free equilibrium and cause treatment failure. However, due to the lack of differentiability at the tumour-free equilibrium, the numerical simulation conducted by Wei [42] did not provide the result of the stability property of the tumour-free fixed point. Furthermore, it is not certain whether the numerical results obtained near the tumour-free fixed point are intrinsic or produced by numerical error.
Evidence suggests that many of us have microscopic tumour cells that remain asymptomatic and never develop or progress to become large and lethal [3,8,20,21,30,32,33,37]. At this stage, a balance between the immune system and the tumour cells could result in long-term persistence of tumour dormancy [3,20,21]. A mathematical model that describes this phenomenon shall have a stable equilibrium of microscopic tumour size. The ratio-dependent functional response assumes that the search time is neglectable. Cytotoxic T cells are able to destroy tumour cells as soon as the tumour cells appear. In this paper, a mathematical model based on the model by de Pillis et al. [18] is proposed where a Beddington-DeAngelis-like functional response is incorporated. The Beddington-DeAngelis functional response considers the search time, handling time, and interference time among predators. When it takes time for cytotoxic T cells to seek and identify tumour cells, the process of tumour cell lysis is delayed and thus tumour cell survival increases. The model with Beddington-DeAngelis-like functional response is expected to have a stable equilibrium of a microscopic tumour size. Furthermore, this model is differentiable. Because of the differentiability at the tumour-free equilibrium, this model can be used to verify whether or not the phenomena generated by numerical simulation of the previously published model [42] are intrinsic. Furthermore, higher dimensional bifurcation analysis than that in the work by Wei [42] will be performed to obtain a better understanding of the system.
To provide details, the mathematical model and mathematical analysis of the microscopic tumour equilibrium are given in Section 2. Some numerical examples and a discussion are given in Section 3. Section 4 provides a brief conclusion.

Mathematical model
The dimensional form of the mathematical model proposed by de Pillis et al. [18] is as follows:Ṫ The state variables T(t), N(t), L(t), C(t), M(t), and I(t) represent the tumour cell population, NK cell population, CD8 + T cell population, number of circulating lymphocytes, chemotherapy drug concentration in the bloodstream, and immunotherapy drug concentration in the bloodstream, respectively. The physical meanings of the terms in Equations (1)-(7) are adapted from de Pillis and Radunskaya [16] and de Pillis et al. [18]. Tumour growth is modelled with a logistic growth equation. The two terms cN and D in Equation (1) represent the fraction of tumour cells killed by NK and CD8 + T cells, respectively. In Equation (2), the growth term for NK cells is proportional to the population of circulating lymphocytes, and the natural death rate is f. The immune cell recruitment terms are gT 2 N/(h + T 2 ) and jD 2 T 2 L/(k + D 2 T 2 ). The terms pNT and qLT represent the inactivation of immune cells after some amount of interaction with tumour cells. Cell growth for CD8 + T cells consists only of the natural death rate m. The terms r 1 NT and r 2 CT describe that CD8 + T cells are also recruited by the debris from tumour cells killed by NK cells and stimulated by the presence of tumour cells. The term uNL 2 describes that regulation and suppression of CD8 + T cell activity occur when there are very high levels of CD8 + T cells. Circulating lymphocytes are generated at a constant rate α, and the natural death rate is β. The terms v M (t), v I (t), and v L (t) are the administration of chemotherapy drug, IL-2, and CD8 + T cells, respectively. The chemotherapy drug follows an exponential decay, −γ M, after injection. Similarly, the immunotherapy drug, IL-2, decays exponentially with −μ I I.
The function D in Equation (7) depends on the relative populations L/T because the experimental data in [19] show that per cent lysis appears to be a function of the ratio of CD8 + T cells to tumour cells [16]. The term DT in Equation (1) is equivalent to DT = dL l T/(sT l + L l ). When l = 1, DT is the ratio-dependent functional response, which has been used in predator-prey models [1,2,4,5,14]. The ratio-dependent functional response considers the handling time and the time wasted due to predator interference. These effects are interpreted as the time it takes for CD8 + T cells to kill tumour cells and the time wasted due to the encounters between CD8 + T cells. A more general functional response is the Beddington-DeAngelis functional response [7,15] which considers the search time, handling time, and predator interference. When the search time is neglectable, the Beddington-DeAngelis functional response becomes the ratio-dependent functional response. Hence DT assumes that the time it takes for CD8 + T cells to seek or recognize tumour cells is neglectable. In this paper, the search time is incorporated in tumour cell lysis by CD8 + T cells, and Equation (7) is modified into where w is related to the search time.
Note that D 1 = D when w = 0, and the tumour-free equilibrium in Equations (1)-(7) is asymptotically stable [42]. When the search time is positive, the effectiveness of tumour cell lysis decreases. Other parameters involved in Equation (8) are d and s representing the saturation level of fractional tumour cell kill by CD8 + T cells and steepness coefficient of the tumour lysis term D, respectively [18]. The parameter value for d does not change in Equation (8) because both D and D 1 approach d when L approaches infinity. Rewrite Equation (8) as follows: The steepness coefficient of tumour lysis becomes w/T l + s for D 1 . When w/T l s, the values for the steepness coefficient are similar for D and D 1 . For example, if we take the parameter values s = 0.512, and l = 1.81 from [18] and let w = 10 −8 , Figure 1 shows that the curves D and D 1 are very close when T = 10 −4 and their steepness are almost the same. The steepness coefficient is larger and the curve D 1 is less steep when T = 10 −5 . Note that T = 10 −4 corresponds to a tumour with 10 5 tumour cells. Using the convention for the cell number to volume conversion where a tumour reaching the size of 1 cm 3 contains 10 9 cells [13], T = 10 −4 (10 5 tumour cells) represents a tumour of size less than 1 mm in diameter. Thus if we choose a small value for w, we can still use the same parameter values for s, d, and l from the previous study [18].
We shall expect that the system using D 1 has a stable small-tumour equilibrium when w takes a small value. Furthermore, the system becomes differentiable at the tumour-free equilibrium, and the local stability analysis is valid for determining the stability of the tumour-free equilibrium. A previous study of Equations (1)-(7) by Wei [42] found that Equation (6) does not affect the dynamics of the system. Also, Equation (6) describes the administration of Interleukin-2 (IL-2), which is immunotherapy treatment. Thus, the effect of Equation (6) can be incorporated into the immunotherapy treatment term v L (t) in Equation (3) for simplicity. Using the following transformations as suggested by de Pillis et al. [18]: and dropping the stars for notational clarity, the non-dimensional form of the model to be studied in this paper is as follows:

Mathematical analysis
The tumour-free case is the basic configuration of the mathematical model because it determines whether the immune system is able to control the growth of a small tumour without treatment and whether the tumour will relapse after attempts to remove the main tumour mass [35]. In this case, only Equation (9)-(12) with Equation (14) are considered and we set v L (t) = 0 and M(t) = 0 for all t ≥ 0. The tumour-free equilibrium is given by Because C is independent of T, N, and L, it may be assumed that C(t) is kept at its equilibrium value 1/β. Thus, the ODE system considered in this subsection becomes as follows:Ṫ where C = 1/β. It is clear that T(t), N(t), and L(t) are non-negative if T(0), N(0), and L(0) are non-negative.
Proof: From Equation (15), The tumour-free equilibrium of Equations (15)- (17) is denoted by T f = (0, 1/f β, 0), and the linearization of the system at T f is The tumour-free equilibrium is stable if c > f β. This condition is equivalent to cN > 1, where N = 1/fg is the number of NK cells of the tumour-free equilibrium. A small tumour will be eliminated if the fraction of tumour cells killed by NK cells is larger than one. However, the parameter values estimated from experimental data of mice and two patients [18] do not satisfied this condition. When the tumour-free equilibrium is unstable, we are interested in whether or not the system has a stable small-tumour equilibrium. The analytical formula for the positive equilibrium is not yet possible, so we will derive parameter conditions for the possibility that the system has a stable equilibrium of small tumour size.

Theorem 2.2:
Assume that f > g and d> 1. and then either there is t 0 > 0 such that T(t 0 ) < 0 or T(t) ≤¯ for t > 0.
Theorem 2.2 provides a prediction of whether or not an immune system might be able to control a small tumour. From Equation (21), Equations (20) and (23) give a range of possible 0 and¯ values satisfying the conditions in Theorem 2.2. Table 1 shows two sets of parameter values in Equations (9)- (14). These parameter values are provided by de Pillis et al. [18]. Note that the data for the search time are not available. Evidence has shown that microscopic tumours (many of them less than 1 mm in diameter) are found in autopsy on individuals who have not been diagnosed with cancer at the time of death [32]. From Figure 1, if we set w = 10 −8 , the search time affects the efficiency of tumour cell lysis only when T < 10 −4 . This may cause changes in other parameter values because the functional response has been changed in the model. However, using the published parameter values with w = 10 −8 allows us to compare the results from the current model and previously publish model and helps understand the effect of the search time on the control of tumour growth. Parameter values of Patient 9 with 0 = 5 × 10 −4 and = 10 −3 satisfy the conditions in Theorem 2.2. It is expected that the immune system of Patient 9 is able to control a small tumour. Parameter values of Patient 10 with a smaller m value, for example m < 18, and the same 0 and¯ values as above also satisfy the conditions.

Bifurcation analysis for immune system response to tumour
From the analysis result obtained in Section 2.2, Equations (15)-(18) are expected to be able to produce a stable equilibrium representing a microscopic tumour, which has a size less than 2 mm in diameter [37]. To study the dynamics of Equations (15)-(18), we begin with a three-dimensional bifurcation diagram as shown in Figure 2(a). Recently, numerical methods for bifurcation analysis have been constructed for ODE systems [39,40], and these numerical methods are applied to Equations (15)- (18). Based on Equations (1)-(7), both NK cells and CD8 + T cells are capable of killing tumour cells. However, the effectiveness c, shown in Table 1, of NK cells in lysing tumour cells is so small that the term cNT does not affect the dynamics of the system. Whether or not the immune system is able to control a small tumour mainly relies on the lysis by CD8 + T cells.
In Figure 2(a), d, m, and β are used as bifurcation parameters, which represent the effectiveness of tumour cell lysis by CD8 + T cells, the death rate of CD8 + T cells, and the death rate of circulating lymphocytes, respectively. Other parameter values are taken from patient 10 in Table 1. When β = 0, the number of circulating lymphocytes grows unboundedly. The case β = 0 is excluded in the numerical simulation. The computational domain is (m, d, β) ∈ [0, 25] × [0, 7] × [0.001, 1]. The system has 4 possible equilibria, labelled by E 0 (the tumour-free equilibrium), E 1 (the stable small-tumour equilibrium), E 2 (the unstable tumour equilibrium), and E 3 (the stable large-tumour equilibrium). Each region has a different asymptotic behaviour. The stable equilibria are labelled in each region. The system with parameter values lying under the surface exhibits bistability, where E 1 and E 3 are stable and E 0 and E 2 are unstable. The immune system is able to control a tumour with a restricted initial size. As the system moves above the surface, it undergoes a fold bifurcation where E 1 and E 2 coincide and disappear. The equilibrium E 3 becomes globally stable. This situation occurs when the death rate of CD8 + T cells or circulating lymphocytes increases. Recall that the term r 2 CT represents the recruitment of CD8 + T cells at the presence of a tumour. When the death rate of circulating lymphocytes is large, the equilibrium number C of circulating lymphocytes is small. The immune system is not able to recruit enough CD8 + T cells to fight tumour cells. A similar result occurs when the death rate of CD8 + T cells is large. In this case, CD8 + T cells die out quickly. The birth rate of tumour cell population in the non-dimensional model, Equations (15)- (18), is equal to one. If the maximum level of fractional tumour cell kill d < 1, a tumour grows to almost its carrying capacity eventually. Theorem 2.2 shows that d > 1 is a necessary condition for the immune system to be able to control a small tumour. Figure 2(a) confirms this condition.
A two-dimensional bifurcation diagram using m and β as bifurcation parameters is shown in Figure 2(b). In this case, we let d = 4.36 where CD8 + T cells are effective in lysing tumour cells. Figure 2(b) shows that the fold bifurcation curve is located near the positive axes indicating that at least one of the death rates, m and β, has to be small for the immune system to have enough CD8 + T cells to control a small tumour. Whether or not the immune system is able to control a small tumour depends on individual variation. Figure 2(c) shows a phase-parameter bifurcation diagram using m as the bifurcation parameter and we let β = 0.0186. The T axis has been rescaled for clarity purposes. Theorem 2.2 predicts that the immune system is able to control a small tumour if the death rate of CD8 + T cells satisfies m < 18. Figure 2(c) gives a similar but larger range (m < 23.6). The small-tumour equilibrium represents a microscopic tumour of 6.3 × 10 3 to 3 × 10 5 cells, where T = 1 is equivalent to 10 9 cells. The unstable equilibrium curve for E 2 can be used to predict the maximum tumour burden, ranging from 3 × 10 5 to 8 × 10 6 cells, which the immune system is able to control with respect to the death rate m of CD8 + T cells. As the death rate m increases, the maximum tumour burden E 2 that the immune system is able to control decreases and the equilibrium tumour size E 1 increases. A fold bifurcation occurs at m = 23.6 where E 1 and E 2 coincide and disappear. This m value gives a death rate threshold for CD8 + T cells. The large-tumour equilibrium E 3 becomes globally stable when m > 23. 6 indicating that a small tumour grows to its carrying capacity eventually. Numerical simulation using parameter values for patient 9 gives a bifurcation diagram, which is not shown here, similar to Figure 2(c) indicating that the immune system of patient 9 is able to control a microscopic tumour.

Bifurcation analysis for chemotherapy treatments
For the numerical simulation of cancer treatment, the treatment schedule in the work by Wei [42] is adopted. In this treatment schedule, chemotherapy starts a week before immunotherapy, and it is followed by an intravenous infusion of TIL on day 8. The cycle of treatment is τ = 9.051, which corresponds to 3 weeks. The forms of v M and v L are summarized as follows: To simulate the outcome of a treatment using chemotherapy alone, we let d L = 0. The parameter values of patient 10 are used in this numerical simulation. The ODE system, Equations (9)- (14), is turned into a discrete dynamical system using a τ -shift map. Then, the numerical method for bifurcation analysis constructed by Wei [41] is directly applied to the τ -shift map. Figure 3  stable for m < 23.6. A small dosage of a chemotherapy drug takes the system to the region where E 3 is globally stable unless the death rate of CD8 + T cells is small. A large dosage is required to eliminate a tumour, and the system is in the region where E 0 is globally stable. As the dosage d M increases, the system undergoes three fold bifurcations and one transcritical bifurcation when the death rate m ∈ (1.5, 18.5).
To better understand these bifurcations and the changes in the equilibria, we let the death rate m = 3. Figure 3(b) shows the bifurcation diagram using the dosage d M as the bifurcation parameter. When d M = 0, both E 1 and E 3 are stable, and the immune system is able to control a tumour with T < 0.006. As the dosage d M increases, the equilibrium tumour size increases. The system undergoes a fold bifurcation at d M = 0.625, where E 1 and E 2 coincide and disappear. As the dosage d M increases past 0.625, the large-tumour equilibrium E 3 becomes globally stable, and the immune system fails to be able to control a small tumour, which would remain microscopic if untreated. The chemotherapy drug deteriorates the host's immune system. This also shows that chemotherapy-induced damage to the immune system observed in the previous model [42] is indeed generated by the mathematical model. As the dosage d M continues to increase, the system undergoes the second fold bifurcation at d M = 9.4, where E 1 and E 2 appear. The chemotherapy regimen is able to treat a tumour smaller than E 2 at a dosage d M > 9.4. As d M increases further, the equilibrium population size E 1 decreases, and the system undergoes the third fold bifurcation at d M = 18.2, where E 2 and E 3 coincide and disappear. The small-tumour equilibrium E 1 becomes globally stable, and the treatment is able to reduce the size of a large tumour to a microscopic size. The system undergoes a transcritical bifurcation at d M = 21, where E 0 and E 1 coincide. The tumour-free fixed point E 0 becomes globally stable if d M > 21, and the treatment is able to eliminate a tumour of any size. Note that the chemotherapy strength ranges from 2 to 5 in the numerical simulation by de Pillis et al. [18]. This range is equivalent to a dosage range from 4.6 to 11.6 for the non-dimensional parameter d M .
Therefore, the above simulation shows that it requires a large dosage for a chemotherapy treatment to be effective.

Bifurcation analysis for combination treatments
In this subsection, we study the effect of combination treatments using immunotherapy and chemotherapy. Figure 4(a) shows a bifurcation diagram using d L and d M as bifurcation parameters, which represent the dosages of immunotherapy drug and chemotherapy drug, respectively. The computational domain is 25]. Figure 4(a) contains a fold bifurcation curve and a transcritical bifurcation curve. When the dosage d M = 0, increases in the dosage d L do not affect the long-term tumour state, meaning that immunotherapy alone is ineffective. The small-tumour fixed points E 1 and large-tumour fixed points E 3 coexist when d M = 0. As the dosage d M increases, the system undergoes a transcritical bifurcation at d M = 1.1, and the tumour-free fixed point E 0 becomes stable. The mixed therapy is able to eliminate a tumour with a restricted initial size. As the dosage d M increases further, the system undergoes a fold bifurcation at d M ≈ 3.5 for d L > 5, and the tumour-free fixed point E 0 becomes globally stable. A treatment with a dosage d M > 3.5 is able to eliminate a large tumour. In contrast, Figure 3(a) shows that a treatment using chemotherapy alone requires a large dosage d M > 21 for E 0 to be globally stable. Increases in the dosage d M do not produce adverse effect. Combination treatments reduce or prevent chemotherapy-induced damage which occurs only when the dosage d L is very small. Figure 4(b). shows the magnification of the portion 0 ≤ d L ≤ 1.5 × 10 −3 in Figure 4(a). Chemotherapy-induced damage does not occur if d L > 6 × 10 −4 , and so this damage can be prevented with a small number of TIL cells.
In the next example, we let d L = 2, which corresponds to about 10 9 TILs.

Conclusion
This paper studies the dynamics of a mathematical model, based on a previously published model developed by de Pillis et al. [18], of tumour growth with mixed immunotherapy and chemotherapy of cancer. The previously published model assumes that the search time is neglectable meaning that CD8 + T cells are able to attack tumour cells immediately after tumour cells appear. The previously published model contains a stable tumour-free equilibrium. The immune system is able to eliminate a small tumour depending on the strength of the immune system. In this paper, a positive search time is incorporated into the model assuming that it takes time for the immune system to seek or recognize tumour cells before the immune system starts to destroy them. The tumour-free equilibrium of the system is unstable, and a stable microscopic tumour equilibrium, called cancer without disease, exists when the search time is small.
An estimate of the stability condition is derived for the microscopic-tumour equilibrium. We have found that the parameter values for patient 9 satisfy the stability condition, and the parameter values for patient 10 satisfy the stability condition with a slightly smaller death rate for CD8 + T cells. Numerical simulation has shown that the microscopic-tumour equilibrium is stable for both patient 9 and 10 meaning that both patients are able to control a microscopic tumour without treatment (see Figure 2(c)). From the mathematical model, the death rate of circulating lymphocytes and the death rate of CD8 + T cells significantly influence the stability of the microscopic-tumour equilibrium because they significantly influence the number of CD8 + T cells for fighting cancer. The immune system is more likely to be able to control a microscopic tumour if the host has a small death rate for circulating lymphocytes or CD8 + T cells. An increase in either one of these two death rates causes the large-tumour equilibrium to become globally stable (see Figure 2(b)). In this case, a tumour will grow to its carrying capacity if not treated.
It is well known that chemotherapy destroys not only the cancer cells but also the immune system responses. The numerical result obtained in Figure 3(b) shows that the use of chemotherapy drugs may damage the immune system causing a microscopic tumour, which would otherwise remain microscopic if not treated, to grow large. A large dosage is required for the tumour-free equilibrium to become locally stable, and a much larger dosage is needed for it to become globally stable (see Figure 3(a)). The numerical simulation, Figure 4(a), also shows that immunotherapy alone is ineffective; increases in TIL cells do not affect the dynamics of the system. A combination of immunotherapy and chemotherapy is then considered for the treatment terms in Equations (9)- (14). The combination therapy using immunotherapy and chemotherapy produces a synergistic effect and allows the use of a small dosage of each therapy drug for successful treatment (see Figure 4(c)). Furthermore, the combination therapy can prevent the chemotherapyinduced damage (see Figure 4(b)). This agrees with clinical studies which have suggested that the combination therapy using chemotherapy and immunotherapy may overcome the chemotherapy-induced lymphocyte damage and improve the outcomes of treatments [6,23,28]. Further studies including clinical tests are needed to verify these results.