Abstract
When it comes to improving cancer therapies, one challenge is to identify key biological parameters that prevent immune escape and maintain an equilibrium state characterized by a stable subclinical tumor mass, controlled by the immune cells. Based on a space and size structured partial differential equation model, we developed numerical methods that allow us to predict the shape of the equilibrium at low cost, without running simulations of the initial-boundary value problem. In turn, the computation of the equilibrium state allowed us to apply global sensitivity analysis methods that assess which and how parameters influence the residual tumor mass. This analysis reveals that the elimination rate of tumor cells by immune cells far exceeds the influence of the other parameters on the equilibrium size of the tumor. Moreover, combining parameters that sustain and strengthen the antitumor immune response also proves more efficient at maintaining the tumor in a long-lasting equilibrium state. Applied to the biological parameters that define each type of cancer, such numerical investigations can provide hints for the design and optimization of cancer treatments.
Graphical Abstract
1 Introduction
The immune system plays a major role in the control of tumor growth. This has led to the concept of immune surveillance and cancer immunoediting composed of three phases (–): the elimination, when tumors are rapidly eradicated by the immune system, the equilibrium, a latency period when tumors can survive but remain on a controlled state, and the escape, the final outgrowth of tumors that have outstripped immunological restraints. In this later phase, immune suppression is prevailing and immune cells are also subverted to promote tumor growth. Numerous cancer immunotherapy strategies have been designed and assessed to counteract immune suppression and restore effective and durable elimination of tumors (–). They show improved efficacy over conventional anticancer treatments but only a minority of patients respond. The challenge to face now is to identify key biological parameters which will convert a fatal outcome into a chronic, manageable state, the durable maintenance of cancer in a viable equilibrium phase controlled by immunity. Reaching such immune-mediated tumor mass dormancy is indeed the first key step for successful control of tumor growth and a goal for immunotherapy (). The equilibrium state is however difficult to apprehend experimentally because the tumor mass at equilibrium is below detectable limits (). Mathematical modeling of the tumor-immune system interactions offers useful information about the features of the equilibrium phase during primary tumor development, and such tools could be used to guide the design of optimal anticancer therapies (–).
We previously () introduced a specific multiscale mathematical model based on partial differential equations (PDE), intended to describe the earliest stages of tumor-immune system interactions. We conjecture that the space heterogeneities of the distribution of active and resting immune cells, which are subjected to several interaction mechanisms with the tumor cells, plays a critical role in the efficiency of the immune response, and the ability in reaching the equilibrium phase. This, in turn, motivates the appeal to PDEs descriptions and can complete the already established modeling based on ordinary differential systems, on which there exists a wide literature, see for instance (, –) Extension to the PDE framework has permitted to bring out the role of space organisation (–). The reader can find further details and references about the mathematical modeling of tumor-immune system interactions, based on different viewpoints and addressing several issues of the efficacy of the immune response, in the reviews (–). The original model developed in () thus accounts for both the growth of the tumor, by natural cell growth and cell divisions, and the displacement of the immune cells towards the tumor, by means of activation processes and chemotaxis effects. The most notable finding from () was that an equilibrium state, with residual tumor and active immune cells, can be observed. Moreover, mathematical analysis provides a basis for the explanation of the formation of the equilibrium. How the biological parameters shape this equilibrium is the main question investigated in the present article. Indeed, the equilibrium can be mathematically interpreted by means of an eigenproblem coupled to a stationary diffusion equation with constraint. This observation permits us to develop an efficient numerical strategy to determine a priori the shape of the equilibrium — namely, the size distribution of the tumor cells and the residual tumor mass — for a given set of biological tumor and immune cell parameters. Consequently, the equilibrium state can be computed at low numerical cost since we can avoid the resolution of the evolution problem on a long time range. The use of this simple and fast algorithm allows us to address the question of the sensitivity of the residual mass to the parameters and to discuss the impact of treatments. This information can be decisive to design clinical studies and choose therapeutic strategies that will revert to an equilibrium phase. Our work therefore provides hints for cancer treatment management.
1.1 Quick Guide to Equations: A Coupled PDE Model for Tumor-Immune System Interactions
The modeling approach imposes to select a few phenomena, considered as the leading effects for the situation under consideration; other effects are just roughly described by tuning some parameters or are simply disregarded. Choices for designing the mathematical model are also dictated by the difficulty in attributing numerical values to the parameters of the equations, due to a lack of experimental measurements: the poor knowledge of driving quantities leads to keep a description as simple as possible, with a reduced number of unknown parameters. The principles of the modeling adopted in (), summarized by Figure 1, led to couple an evolution equation for the size-distribution of the tumor cells, and a convection-diffusion equation for the activated immune cells. The two-way coupling arises from the death term induced by the action of the immune cells on the tumor cells, and by the activation and the attraction of immune cells towards the tumor, which are determined by the total mass of the tumor. The model is intended to describe the earliest stages of the tumor formation, when the size of the tumor is relatively small. The tumor is located at the center of a domain Ω (there is no displacement of the tumor). The model distinguishes two distinct and independent length scales: the size of the tumor cells, described by the variable z≥0 , is considered as “infinitely small” compared to the scale of displacement of the immune cells, described by the space variable x ∈ Ω .
Figure 1
The unknowns are
The size density of tumor cells ( t, z )↦n( t, z ) so that the integral gives the volume of the tumor occupied at time t by cells having their size z in the interval (a, b);
The concentration of activated immune cells which are fighting against the tumor ( t, x )↦c( t, x ) ;
The concentration of chemical signal that attracts the immune cells towards the tumor microenvironment ( t, x )↦ϕ( t, x ).
The specific biological assumptions made to construct the model are fully described in (). Figure 2 offers an overview of the interaction mechanisms embodied in the equations and of the role of the parameters of the model.
Figure 2
Immune cells, once activated from a bath of resting cells, are subjected to natural diffusion and to a chemotactic drift, induced by the presence of the tumor. The strength of this drift, as well as the activation of immune cells, directly depends on the total mass of the tumor, proportional to the quantity
The immune system-tumor competition is described by the following system of PDEs
The features of the growth-division dynamics for the tumor cells (1a) are embodied into the (possibly size-dependent) growth rate z↦V( z )≥0 and the cell division operator Q(n). We refer the reader to (–) for further details on this evolution equation (with m (n, c) =0) for cell growth and division, and its application to cancer modeling. What is crucial for modeling purposes is the principle that cell-division does not change the total mass: the operator Q satisfies . However, the total number of cells in the tumor increases since (we refer the reader to () and Appendix A for further details). In what follows, we restrict to the mere symmetric binary division operator
with z↦a(z )≥0 the division rate. It simply describes the situation where cells are cut into two cells having half the size of the original cell. Further relevant examples of division operators can be found in () (see Appendix A). The specific case where the division rate a in () is a positive constant makes the model simpler, and is often used. It is however likely relevant to incorporate more complex behaviors through the size-dependence; for instance divisions can be prohibited below a certain size threshold. Similarly, it can be convenient to assume that the growth rate V is a positive constant, but more intricate laws can take into account some important phenomena. For instance, logistic or Gompertz law can incorporate size limitation effects, and roughly describe difficulties in accessing nutrients or necrotic effects (–); a detailed study of growth laws can be found in (). As mentioned above, though, using such complex laws, also raises the issue of determining more parameters. The boundary condition for n in (1d) means that no tumor cells are created with size 0.
Despite the fact that there exists several types of immune cells – at least T-cells and NK cells – fighting against the tumor, they are all described here through the single concentration c. It also means that coefficients of the equation – the death rate γ>0, the chemotactic strength χ>0 , and the diffusion coefficient D – correspond to an averaged behavior of all these cells. By the way, working with a constant diffusion coefficient D > 0 is again a simplification, neglecting the architecture of the tumor environment, which might induce directional effects. The effector immune cells that effectively fight against the tumor, are activated from a “reservoir” of resting cells, described in the right hand side of (1b) by ( t,x )↦R( t,x ) . This given function, possibly time and space dependent, stands for the space distribution of the influx rate of activated effector immune cells. It takes into account the sources of resting immune cells that can be activated in the tumor microenvironment or in the draining lymph nodes into cells fighting the tumor. At early stages of tumor growth, the rate of the activation process is supposed to be directly proportional to the tumor mass μ1. Again, more complex activation law, for instance based on Michaelis-Menten kinetics can incorporate relevant limitation mechanisms. The Dirichlet boundary condition for c in (1d) means that the immune cells far from the tumor are non-activated. Immune cells are directed towards the tumor by a chemo-attractive potential ϕ, induced by the presence of the tumor cells. Through (1c), the strength of the signal is proportional to the total mass of the tumor, and it is shaped by a form function x↦σ( x ) which will be a function peaked at the tumor location. The potential is thus defined by the diffusion equation (1c), that involves a positive coefficient K>0 (that could be matrix valued), and the Neumann boundary condition in (1d), where v stands for the unit outward normal vector on ∂Ω. Finally, the activated immune cells are able to destroy tumor cells, as described by the death term in (1a)
where δ≥0 is another form function, also peaked in the vicinity of the tumor. For the numerical experiments, we shall work with the Gaussian profiles
where the positive parameters A,Aσ and θ,θσ can be used to tune the amplitude and spreading of these functions, and thus the strength and radius of influence of the related phenomena. We refer the reader to () for further details and comments about the model. Note that this model neglects the possible additional protumoral effects that can take place and are crucial to swing to the escape phase. Such protumor effects can have different forms: they can directly enhance the tumor growth, and make antitumor immune cells exhausted, a state where they are hyporesponsive and cannot kill the tumor, see () on these issues. Remarkably, the model (1a)-(1e) is able to reproduce equilibrium phases where the tumor growth is controlled by the immune response.
2 Materials and Methods
2.1 Development of Numerical Methods Predicting Parameters of the Equilibrium in Immune-Controlled Tumors
According to (, , ), the equilibrium phase corresponds to a long-lasting period of immune-mediated latency, also known as tumor mass dormancy, prior to the emergence of clinically detectable malignant disease, with a residual tumor which has not be fully destroyed by the immune system, maintained under the control of immunity. The simulations of the initial-boundary value problem (1a)-(1e) performed in () revealed that such a behavior can be reproduced by the model. Here, we wish to study the features of the equilibrium phase in immune-controlled tumors and, in particular, we want to predict, for given biological parameters (see Section 2.2 below), the total mass of the residual tumor and its size distribution. To this end, we developed specific numerical procedures based on the mathematical interpretation of the equilibrium.
2.1.1 Equilibrium States
The definition of the equilibrium relies on the following arguments. When disregarding the immune response, the cell-division equation
admits a positive eigenstate, which drives the large time behavior of the solution. To be more specific, there exists λ>0 and a non negative function satisfying
The existence-uniqueness of the eigenpair can be found in (, ). Furthermore, when the tumor does not interact with the immune system, the large time behavior is precisely driven by the eigenpair: the solution of () behaves like
where μ0 >0 is a constant determined by the initial condition, see (, ). Consequently, in the immune-free case, the tumor population grows exponentially fast, with a rate λ>0 , and, as time becomes large, its size repartition obeys a certain profile . In the specific case where V is constant and Q is the binary division operator (), with a constant division rate a, we simply have λ=a and the profile is explicitly known (, ). However, for general growth rates and division kernels the solution should be determined by numerical approximations; we are going to detail a numerical procedure to effectively compute the pair .
Coming back to the coupled model (1a)-(1e), we infer that the equilibrium phase corresponds to the situation where the death rate – the integral of the immune cells concentration with weight δ, denoted as in (3) – precisely counterbalances the natural exponential growth of the tumor cell population. In other words, at equilibrium we expect that
• The size distribution of tumor cells is proportional to the eigenstate . The proportionnality factor is related to the total mass by the relation .
• The concentration of immune cells is defined by the stationary equation
• Where Ф is the solution of
• Endowed with the homogeneous Neumann boundary condition, together with the constraint
This can be interpreted as an implicit definition of the total mass μ1 to be the value such that the solution of the boundary value problem (7) satisfies (8): in other words, it defines implicitly the mass of the residual tumor μ1 to be the value such that the solution of the stationary boundary value problem for C defines a death rate that exactly compensates the exponential growth rate of the growth division equation. The existence of an equilibrium state defined in this way is rigorously justified in (10, Theorem 2).
Theorem 2.1. Let x↦R( x )∈L2 ( Ω ) be a non negative function. If λ>0 is small enough, there exists a unique μ1 ( λ )>0 such that the solution Cμ1 ( λ ) of the stationary equation (7) satisfies (8).
Theorem 2.1 requires a smallness assumption; for (2) with constant growth rate V and division rate a, this is a smallness assumption on a. Numerical experiments have shown different large time behaviors for the initial-boundary value problem (1a)-(1e) (an example will be presented later on):
When the source term R is space-homogeneous, the expected behavior seems to be very robust. The immune cell concentration tends to fulfill the constraint as time becomes large, and the size repartition of tumor cells tends to the eigenfunction . The total mass μ1 tends to a constant; however the asymptotic value cannot be predicted easily.
When R has space variations, the asymptotic behavior seems to be much more sensitive to the parameters of the model, in particular to the aggressiveness of the tumor (characterized by the cell division rate a). On short time scale of simulations, we observe alternance of growth and remission phases, and the damping to the equilibrium could be very slow.
These observations bring out the complementary roles of different type of cytotoxic cells (). The NK cells could be seen as a space-homogenous source of immune cells, immediately available to fight against the tumor, at the early stage of tumor growth. In contrast, T-cells need an efficient priming which occurs in the draining lymph nodes, and their sources is therefore non-homogeneously distributed. Eventually, NK and CD8+ T-cells cooperate to the anti-tumor immune response.
Numerical experiments thus show that the model (1a)–(1e) is able to reproduce, in the long-time range, cancer-persistent equilibrium, but the features of the equilibrium, and its ability to establish, are highly sensitive to the parameters. To discuss this issue further, we focus here on the mass at equilibrium considered as a critical quantity that evaluates the efficacy of the immune response. Indeed, it is known that a tumor gains in malignancy when its mass reaches certain thresholds (, ). The smaller the tumor mass at equilibrium, the better the vital prognosis of the patient. In doing so, we do not consider transient states and time necessary for the equilibrium to establish. The interest of the interpretation of the equilibrium by means of an eigenproblem relies on the fact that the equilibrium state can be determined a priori, at least through numerical simulations, without running the initial boundary value problem over long time ranges: given a set of biological parameters it can be obtained by solving the eigenvalue problem for and the constrained stationary drift-diffusion equation for C, see Figure 3. In turn, since the equilibrium state can be computed at low numerical cost, a wide range of parameters can be considered and the role of the parameters can be investigated in details. The determination, on numerical grounds, of the equilibrium state relies on a two-step process, as schematised in Figure 3. First, we compute the normalized eigenstate of the tumor cell equation, second, we find the tumor mass which makes the coupled death rate fit with the eigenvalue. To this end, we have developed a specific numerical approach.
Figure 3
2.1.2 The Eigen-Elements of the Growth-Division Equation
The numerical procedure is fully detailed and analyzed in Appendix A; it is inspired from the spectral analysis of the equation: λ is found as the leading eigenvalue of a conveniently shifted version of the growth-division operator. In practice, we work with a problem where the size variable is both truncated and discretized. Hence, the problem recasts as finding the leading eigenvalue of a shifted version of the underlying matrix, which can be addressed by using the inverse power method [(), Section 1.2.5]. We refer the reader to (, ) for a thorough analysis of the approximation of eigenproblems for differential and integral operators, which provides a rigorous basis to this approach. It is also important to check a priori, based on the analysis of the equation (), how large the shift should be, and that it remains independent on the numerical parameters. As already mentioned, for some specific division and growth rates, the eigenpair is explicitly known, see (). We used these formula to validate the ability of the algorithm to find the expected values and profiles.
2.1.3 Computation of the Equilibrium Mass
Having at hand the eigenvalue λ, we go back to the convection-diffusion equation (7) and the constraint (8) that determine implicitly the total mass μ1 of the residual tumor. For a given value of μ1, we numerically solve (7) by using a finite volume scheme, see (, Appendix C). Then, we use the dichotomy algorithm to fit the constraint:
The chemo-attractive potential Ф is computed once for all.
Pick two reference values 0 < μa < μb; the mass we are searching for is expected to belong to the interval (μa, μb).
Set and compute the associated solution Cμ 1 of (7) (the subscript emphasizes the dependence with respect to μ1). Evaluate the discrete version of the quantity I=∫δCμ1 dx−λ
If I < 0, then replace μa by μ1, otherwise replace μb by μ1.
We stop the algorithm when the relative error is small enough.
It is also possible to design an algorithm based on the Newton method. However, this approach is much more numerically demanding (it requires to solve more convection-diffusion equations) and does not provide better results.
2.2 Identification of Biological Parameters
In order to go beyond the qualitative discussion of (10), the model should be challenged with biological data. The PDE system (1a)-(1e) is governed by the set of parameters collected in Table 1. Most parameter values were retrieved from previously published experimental results and we propose an estimation of the remaining parameters R, a, V based on the experimental study performed in () where the development of chemically-induced cutaneous squamous cell carcinoma (cSCC) is investigated.
Table 1
| Symbol | Description | Value and unit | References |
|---|---|---|---|
| X | chemotactic coefficient | 8.64×101−8.64×106mm2·mmol−1·day−1 | (Macrophages) () |
| D | natural space diffusion coef. of the cytotoxic effector cells population | 8.64×10−5−10−3mm2·day−1 | (CD8+ T-cells) (, ) |
| R | the normal rate of influx of effector immune cells | estimated | |
| γ | natural death rate of the tumor antigen-specific cytotoxic effector cells | 2×10−2−1 day−1 | (, , , ) |
| A | strength of the immune response | (–) | |
| K | diffusion coefficient for the attractive potential ϕ | 10−2−1 mm2·day−1 | (, ) |
| Aσ | strength of the chemical signal induced by each tumor cell | 5·10−17−0.625×10−16mmol·−1μm3·day−1 | () |
| a | division rate of the tumor cells | log(a)∼N(log(0.12),0.2)(a in day−1) | estimated |
| V | growth rate of the tumor cells | log(V)∼N(log(816.33),0.51)(V in μm3·day−1) | estimated |
Key model parameters and their biophysical meaning.
Calibrating the parameters of the equations is an issue due to the lack of direct measurements, and the fact that experimental data are obtained at the price of the sacrifice of mice. Consequently, beyond the cost of the experiments, it also means that a time evolution of the quantities of interest is usually not affordable. Therefore, a specific procedure should be developed in order to estimate the parameters from the experimental data points. Since the informations on the parameters are quite poor, we restrict to the case where the coefficients a, V, R are constant, which is also a reasonable assumption when dealing with the earliest stages of the tumor development. In order to identify the parameters, we shall use a degraded version of the equations.
Neglecting the immune response, the tumor growth is driven by (5). As explained above, this leads to an exponential growth of the tumor mass, see (–, ). Let , the total number of tumor cells, and
Integrating (9) with respect to size variable, with integration by parts, and bearing in mind that the cell division operator is mass preserving, we thus get
Next, assuming space homogeneity of the immune cells concentration and neglecting the displacement and the natural death rate of the immune cells, the immune cells concentration is driven by
Based on this simplified dynamics, reduced to (9)-(10), we used the Nonlinear Mixed Effects Modeling (NMEM) in order to estimate the parameters a, V, R from the experimental data. Let N denote the number of mice within the population and the vector of longitudinal measurements for the ith mouse: is a typical observation of the mouse i for a given measurement type k ∈ { 0,1,2} (with (0, 1, 2) referring to ( μ0 ,μ1 ,c ) respectively) at time for i ∈ { 1, … ,N } and i ∈ { 1, … ,N } . We suppose that the statistics of the measurements obeys, for
where is the evaluation of the model at time is the vector of the parameters describing the individual i and the residual error model. The inter-individual variability is described by the combination of fixed effects , which, by definition, are constant within the population and along time, and random effects which explain the inter-individual variability among the mice. The positivity of the parameters is ensured by assuming that the individual parameters follow a log-normal distribution. In other words, the random effects are normally distributed with mean zero and a variance-covariance matrix W. For instance W=diag( ω0 ,ω1 ,ω2 ) where the ωk’s stand for the variance of the parameters a, V, R. Therefore, we have
for k∈{ 0,1,2 } . The error model is assumed to be proportional to the model evaluation and is defined as follows:
Where ϵij ∼N( 0,1 ) represents the statistical model residual errors and b(k) is the proportionality factor measuring the relative amplitude of the errors.
2.2.1 Estimation of the Model Parameters
According to the experimental procedure in (), 5× 105 mSCC38 were injected to each mouse at time t0 = 0. Therefore we fixed the initial number of tumor cells to μ0 ( 0 )=5× 105cells. Assuming that each tumor cell is spherically shaped with a radius 15 μ m, we set μ1 ( 0 )=7.1mm3. The initial concentration of immune cells is fixed to c0 = 0: we suppose that initially there is no effector immune cells (or at least it means that the initial concentration of activated immune cells is negligible compared to the concentration of resting cells). Some data points were censored due to the sacrifice of the individual for flow cytometry cell counting. The censored data points have been handled by Limit Of Quantification (LOQ) censoring (). Let be the finite or infinite censoring interval for mouse i, measurement k and time and
where is the conditional distribution of given . Let us collect in a vector α=( apop,Vpop, Rpop, ωa, ωV, ωR, ba, bV, bR) the parameters of the model; they are estimated by maximizing the observed likelihood function
To this end, we used the Stochastic Approximation of the Expectation Maximization algorithm (SAEM) implemented in the MONOLIX R API (). Furthermore, the individual parameter estimators are computed in MONOLIX () by means of the Empirical Bayes Estimate (EBE) of which corresponds to the mode of the conditional distribution (where corresponds to estimated parameters).
A preliminary estimation procedure indicates a significant correlation between the parameters a and R ( t-test p-value 2.6× 10−6) . Hence, introducing this correlation into the variance covariance matrix of the random effects by setting covar ( a,R )=ρaRωaωR , where ρaR represents the correlation coefficient between a and R, enhances the goodness of fit. The estimated value of ρaR is 0.8 with a relative standard error of 13%. The parameters in a were estimated with reasonable standard errors (computed using the stochastic approximation) and relative standard errors (max (R.S.E) = 30.6 and min (R.S.E.) = 3) which indicate that the model parameters are identifiable. The ShapiroWilk test reinforces the normality hypotheses on the random effects (the p-values for ηa ,ηV and ηR are respectively 0.83, 0.61, 0.2). Pictures indicating the fits are provided in Figure 4, and detailed parameter estimates are given in Table 2.
Figure 4
Table 2
| parameters | value | S.E | R.S.E (%) |
|---|---|---|---|
| apop | 0.12 | 0.0041 | 3 |
| Vpop | 816.33 | 92.59 | 11 |
| Rpop | 2.2×10−6 | 3.6×10−7 | 16 |
| ωa | 0.20 | 0.027 | 13.5 |
| ωV | 0.51 | 0.075 | 15 |
| ωR | 0.84 | 0.11 | 13 |
| ba | 0.37 | 0.041 | 11 |
| bV | 0.17 | 0.052 | 31 |
| bR | 0.18 | 0.056 | 30 |
| ρaR | 0.8 | 0.1 | 13 |
Estimated value of the parameters with their Standard Error (S.E.) and Relative Standard Error (R.S.E).
2.3 Materials
2.3.1 Mice
FVB/N wild-type (WT) mice (Charles River Laboratories, St Germain Nuelles, France) were bred and housed in specific-pathogen-free conditions. Experiments were performed using 6-7 week-old female FVB/N, in compliance with institutional guidelines and have been approved by the regional committee for animal experimentation (reference MESR 2016112515599520; CIEPAL, Nice Côte d’Azur, France).
2.3.2 In Vivo Tumor Growth
mSCC38 tumor cell line was established from DMBA/PMA induced sSCCs and maintained in DMEM (Gibco-ThermoFisher Scientific, Courtaboeuf, France) supplemented with 10% heat-inactivated fetal bovine serum (FBS) (GE Healthcare, Chicago, Illinois, USA) penicillin (100 U/ml) and streptomycin (100 μg/ml) (Gibco-ThermoFisher Scientific, Courtaboeuf, France). 5× 105 mSCC38 were intradermally injected in anesthetized mice after dorsal skin shaving. Tumor volume was measured manually using a ruler and calculated according to the ellipsoid formula: Volume=Length(mm )×Width( mm )×Height( mm )×π/6.
2.3.3 Tissue Preparation and Cell Count
mSCC38 were excised and enzymatically treated twice with collagenase IV (1mg/ml) (Sigma-Aldrich, St Quentin Fallavier, France), and DNase I (0.2 mg/ml) (Roche Diagnostic, Meylan, France) for 20 minutes at 37°C . Total cell count was obtained on a Casy cell counter (Ovni Life Science, Bremen, Germany). Immune cell count was determined from flow cytometry analysis. Briefly, cell suspensions were incubated with anti-CD16/32 (2.4G2) to block Fc receptors and stained with anti-CD45 (30-F11)-BV510 antibody and the 7-Aminoactinomycin D (7-AAD) to identify live immune cells (BD Biosciences, Le Pont de Claix, France). Samples were acquired on a BD LSR Fortessa and analyzed with DIVA V8 and FlowJo V10 software (BD Biosciences, Le Pont de Claix, France).
2.3.4 Mathematical and Statistical Analysis
Computations were realized in Python and we made use of dedicated libraries, in particular the gmsh library for the computational domain mesh generation, the packages optimize (for the optimization methods using the Levenberg-Marquard mean square algorithm; similar results have been obtained with the CMA-ES algorithm of the library cma) from the library scipy, the MONOLIX R API and application for the model calibration to the experimental data (), the library Pygpc for the generalized Polynomial Chaos approximation () and the library Salib for the sensitivity analysis ().
3 Results
3.1 Validation of the Method
For all the simulations discussed here, we adopt the same framework as in (10): the tumor is located at the origin of the computational domain Ω, which is the two-dimensional unit disk. Otherwise explicitly stated, we work with the lower bound of the parameters collected in Table 1. When necessary, the initial values of the unknowns are respectively μ0 (0) = 1 celln, ue μ1 (0) = 14137.2 μm3, c(0,x) = 0.
To start with, we perform a simulation of the initial-boundary value problem (1a)-(1e). Figure 5 illustrates how the equilibrium establishes in time: as time becomes large, the effective concentration of active immune cells, that is denoted
Figure 5
tends to the eigenvalue of the cell-division equation, the total mass μ1(t) tends to a constant and the size distribution of tumor cells takes the profile of the corresponding eigenstate. This result has been obtained by setting ( a,V,R )=( 0.072, 713.61, 1.74× 10−7 ) . We observe a non symmetric shape of the size distribution of tumor cells, peaked about a diameter of 23 μm, which is consistent with observational data reporting the mean size distribution of cancer cells (, ).
For the simplest model of growth-division with a and V constant, we know an expression of the eigenstate ; however, we do not know an explicit evaluation of the residual mass. Nevertheless, we can compare the results of the inverse power-dichotomy procedure that predicts the residual mass, to the large time simulations as performed in (10). Let be the asymptotic value of the total mass given by the large time simulation of the initial-boundary value problem (and checking that the variation of the total mass has become negligible) and let be the mass predicted by the power-dichotomy procedure. We set
The results for several cell division rates a are collected in Table 3: the numerical procedures finds the same equilibrium mass as the resolution of the evolution problem, which is a further validation of the method.
Table 3
| a | (mm3) at final time T = 500 | (mm3) | Eμ 1 |
|---|---|---|---|
| 0.103 | 7.67271875 × 10-5 | 7.67271872 × 10-5 | 4.10 × 10-9 |
| 0.15 | 1.11701535 × 10-4 | 1.11701543 × 10-4 | 7.97 × 10-8 |
| 0.20 | 1.48924575 × 10-4 | 1.48924641 × 10-4 | 4.40 × 10-7 |
| 0.3 | 2.23420663 × 10-4 | 2.23420562 × 10-4 | 4.53 × 10-7 |
| 0.351 | 2.61368442 × 10-4 | 2.61367974 × 10-4 | 1.80 × 10-6 |
Comparison of the large time tumor mass and the predicted tumor mass for several values of a.
Further validation concerning the ability in finding the leading eigenstate are presented in Appendix A. The method has been successfully employed to predict equilibrium state when dealing with complex growth rate and division operator in ().
3.2 Numerical Simulations Show How Parameters Influence Equilibrium
The numerical methods were next used to assess how the parameters influence the equilibrium. In particular, we wish to assess the evolution of the tumor mass at equilibrium according to immune response and tumor growth parameters.
For the numerical simulations presented here, we thus work on the eigenproblem (6) and on the constrained system (7)-(8). Unless precisely stated, the immune response parameters are fixed to the lower bounds in Table 1. The tumor growth parameters are set to = 0.1 day-1, V = 713.61 μm3day-1 and .
The main features of the solutions follow the observations made in (10), which were performed with arbitrarily chosen values for the parameters. We observe that tends to the division rate a, which in this case corresponds to the leading eigenvalue of the cell-division equation. It is important to note that the predicted diameter of the tumor at equilibrium — see Figure 5 — is significantly below modern clinical PET scanners resolution limit, which could detect tumors with a diameter larger than 7mm (). This is consistent with the standard expectations about the equilibrium phase (3), but, of course, it makes difficult further comparison of the prediction with data.
The aggressiveness of the tumor is characterized by the division rate, the variations of which impact the size of the tumor at equilibrium: the larger a, the larger the residual tumor, see Figure 6A. Increasing the immune strength A increases the efficacy of the immune response, reducing the size of the residual tumor see Figure 6B. Similarly, increasing the mean rate of influx of effector immune cells in the tumor microenvironment R, decreases the tumor size at equilibrium, see Figure 6C. On the contrary, increasing the death rate of the immune cells γ reduces the efficacy of the immune response and increases the equilibrium tumor size see Figure 6D.
Figure 6
Moreover, as mentioned above, not only the parameters determine the equilibrium mass, but they also impact how the equilibrium establishes. Figures 7A–C shows what happens by making the tumor cell division rate a vary. There are more oscillations along time, with larger amplitude, as a increases. Similar observations can be made when reducing the strength of the immune system A (likely out of its realistic range), see Figures 7D–F. The smaller A, the weaker the damping of the oscillations and the longer the periods. We notice that the decay of the maximal tumor radius holds at a polynomial rate. In extreme situations, either the damping is very strong and the equilibrium establishes oscillation-free or the equilibrium does not establish on reasonable observation times, and the evolution can be confounded with a periodic alternance of growing and remission phases. Such scenario illustrates that the relevance of the equilibrium can be questionable depending on the value of the parameters. In what follows, we focus on the details of the equilibrium itself, rather than on the transient states.
Figure 7
3.3 Global Sensitivity Analysis on the Equilibrium Mass Identifies the Key Parameters to Target in Cancer Therapy
Since the equilibrium state can be computed for a reduced numerical cost (it takes about 1/4 of a second on a standard laptop), we can perform a large number of simulations, sampling the range of the parameters. This allows us to discuss in further details the influence of the parameters on the residual mass and, by means of a global sensitivity analysis, to make a hierarchy appear according to the influence of the parameters on this criterion. Ultimately, this study can help in proposing treatments that target the most influential parameters.
Details on the applied methods for the sensitivity analysis can be found in Appendix B. Among the parameters, we distinguish:
The tumor cell division rate a which drives the tumor aggressiveness,
The efficacy of the immune system, governed by the mean influx rate of activated effector immune cells R, the strength of the immune response A, the chemotactic sensitivity χ, the death rate γ of the immune cells, and the strength of the chemical signal induced by each tumor cell Aσ
Environmental parameters such as the diffusion coefficients D (for the immune cells) and K (for the chemokine concentration).
We assume that the input parameters except a and R are independent random variables. Due to the lack of knowledge on the specific distribution of most of the parameters, the most suitable probability distribution is the one which maximizes the continuous entropy (), more precisely, the uniform distribution over the ranges defined in Table 1. Therefore, the uncertainty in the parameter values is represented by uniform distributions for the parameters ( A,χ,D,Aσ ,γ,K ) and by log-normal distributions for the parameters a and R. In what follows, the total mass at equilibrium, μ1, given by the power-dichotomy algorithm, is seen as a function of the uncertain parameters:
To measure how the total variance of the output μ1 of the algorithm is influenced by some subsets i1 ⋯ip of the input parameters i1 ⋯ip ( k≥p being the number of uncertain input parameters), we compute the so-called Sobol’s sensitivity indices. The total effect of a specific input parameter i is evaluated by the total sensitivity index , the sum of the sensitivity indices which contain the parameter i. (Details on the computed Sobol indices can be found in Appendix B). The computation of these indices is usually based on a Monte Carlo (MC) method [see (, )] which requires a large number of evaluations of the model due to its slow convergence rate where N is the size of the experimental sample). To reduce the number of model evaluations, we use instead the so-called generalized Polynomial Chaos (gPC) method [see ()]. The backbone of the method is based on building a surrogate of the original model by decomposing the quantity of interest on a basis of orthonormal polynomials depending on the distribution of the uncertain input parameters θ( ω )=( a,A,R,χ,D,Aσ ,γ,K ), where ω represents an element of the set of possible outcomes. Further details on the method can be found in (). For uniform distributions, the most suitable orthonomal polynomial basis is the Legendre polynomials. The analysis of the distribution of μ1 after a suitable sampling of the parameters space indicates that μ1 follows a log-normal distribution. This distribution is not uniquely determined by its moments (the Hamburger moment problem) and consequently cannot be expanded in a gPC [see ()]. Based on this observation, to obtain a better convergence in the mean square sense, we apply the gPC algorithm on the natural logarithm of the output μ1. Typically, ln(μ1) is decomposed as follows:
where ϵ corresponds to the approximation error, and p represents the highest degree of the expansion. Hence, the dimension of the polynomial basis is given by . We reduce the number of model evaluations to 642 runs by constraining also the parameters interaction order to 2. For our purpose, a degree p = 5 gives a better fit (see Figures 8A, B to the original model and the goodness of fit of the gPC algorithm is measured by a Leave One Out Cross Validation (LOOCV) technique (). The resulting LOO error indicates 0.4% prediction error. The Sobol’s sensitivity indices are then computed from the exponential of the surrogate model () by using Monte Carlo simulations combined with a careful space-filling sampling of the parameters space [see (, )]. For the computations, a sample with N=1.8× 106 points has been used in order to get stable second order Sobol indices. Indeed, the sensitivity indices that are needed to discriminate the impact of the input parameters are the first and total Sobol’ sensitivity indices. Here, the analysis revealed a significant difference between some first order Sobol’ indices and their corresponding total Sobol indices, which indicated the importance of computing also the second order Sobol’ indices.
Figure 8
It is important to stress that the obtained results, and the associated conclusions, could be highly dependent on the range of the parameter values. This observation makes the measurement/estimation of the parameters a crucial issue which can be dependent on the type of cancer analyzed.
3.3.1 Efficacy of the Immune Response
The first order Sobol indices represented in Figure 8C indicate that the parameters which impact the most the variability of the immune-controlled tumor mass at equilibrium are:
The strength of the lethal action of the immune cells on the tumor cells A, by far the most influential, and three additional parameters
The influx rate of activated effector immune cells into the tumor microenvironment R.
The natural death rate γ of the effector immune cells,
And the division rate a of the tumor cells.
This result is consistent with the observations made from the numerical experiments above and in (10), showing a prominent role of the immune response which can be enhanced by increasing either A or R, and decreasing γ. That A is the most influential parameter is not that surprising but it is remarkable how far its importance exceeds that of the other parameters. It is also puzzling to see that the chemotactic sensitivity χ, like the strength of the chemical signal induced by each tumor cell Aσ, the space diffusion coefficient of the effector immune cells D and the diffusion coefficient of the chemokines K , have a negligible influence on the immune-controlled tumor mass, see Figure 8C, whether individually or in combination with other parameters. This result is consistent with the necessity for immune cells to be able to effectively kill the tumor cells once they reach the tumor site. The second order Sobol’ indices indicate that the leading interactions are the pairs (A, R), (A,γ), (R, γ), (a,A), (a, R) and (a, γ). Accordingly, in order to enhance the immune response, an efficient strategy can be to act simultaneously on the immune strength A together with the influx rate of activated immune effector cells R. Increasing such influx into the tumor microenvironment by enhancing the activation/recruitment processes leading to the conversion of naive immune cells into activated immune cells potentiate anti-tumor immune responses. Besides, the natural death rate γ of the effector immune cells combined to A and R have an impact, as well as A combined with the division rate of the tumor cells, a.
3.3.2 The Tumor Aggressiveness
The tumor aggressiveness is mainly described by the cell division rate a. The first order Sobol indice indicates that a influences significantly the tumor mass at equilibrium, and we observe that the total Sobol index of a is higher than the individual one. This indicates that this parameter has strong interactions with the others. By taking a look at Figure 8D we remark that a interacts significantly with the parameters A, R, γ. However, the most significant interaction is the one with A. This suggests that combining therapies targeting tumor and immune cells should be more efficient at maintaining immune-mediated tumor mass dormancy ().
3.3.3 Towards Optimized Treatments
Because equilibrium state can be computed for a reduced numerical cost, it allows a large number of simulation to be performed in a minimal time, so that an extensive sampling of the range of the parameters can be tested. The flexibility of the numerical simulations provides valuable tools to assess the efficiency of a variety of therapeutic strategies and select those that sustain a viable equilibrium and prevent cancer relapses after a surgery or a treatment. Figure 9 illustrates how the equilibrium mass is impacted when combining variations of two parameters, namely the immune strength A combined to the tumor cell division rate a, the mean rate of influx of effector immune cells R or the death rate of effector immune cells γ; and the tumor cell division rate a with the death rate γ. Interestingly, a reduction of the tumor mass at equilibrium can be obtained significantly more easily by acting on two parameters than on a single one. For instance, reducing the tumor cell division rate a from 0.35 to 0.1 cannot reduce the diameter of the tumor below.025 mm, with A = 1; while the final size is always smaller when A = 3.95. This observation highlights the interest of combined treatments having such complementary actions. The interest is two-fold: either smaller residual tumors can be obtained by pairing two actions, or the same final tumor size can be obtained with a combined treatment having less toxicity than a mono-therapy.
Figure 9
4 Discussion
Controlling parameters that maintain immune-mediated tumor mass dormancy is key to the successful development of future cancer therapies. To understand how equilibrium establishes and how it is influenced by immune, environmental and tumor-related parameters, we evaluate the tumor mass which tends to a constant at equilibrium. In this study, we make use of the space and size structured mathematical model developed in (10) to provide innovative, efficient methods to predict, at low numerical cost, the residual tumor mass at equilibrium. By means of numerical simulations and global sensitivity analysis, we identify the elimination rate A of tumor cells by immune cells as the most influential factor. Therefore, the most efficient therapeutic strategy is to act primarily on the immune system rather than on the tumor itself. We also demonstrate the need to develop combined cancer treatments, boosting the immune capacity to kill tumor cells (increase A), the conversion into efficient immune cells (increase R), reducing natural death rate of effector immune cells (decrease γ) and reducing the ability of tumor cells to divide (decrease a). The combination of such approaches definitely outperforms the performances of a single action; it permits to maintain the tumor in a long-lasting equilibrium state, far below measurement capabilities.
Generally, therapeutic strategies are designed to target preformed, macroscopic cancers. Indeed, patients are diagnosed once their tumor is established and measurable, thus at the escape phase of the cancer immunoediting process (1). The goal of successful treatments is to revert to the equilibrium phase and ultimately to tumor elimination. Experimental evidence and clinical observations indicate that such equilibrium exists but it is difficult to study and measure, the residual tumor mass being below detection limits (1, 2, 3). It is regarded as “a immune-mediated tumor mass dormancy” when the rate of cancer cell proliferation matches their rate of elimination by immune cells. In human, cancer recurrence after therapy and long periods of remission or detection of low number of tumor cells in remission phases are suggestive of such equilibrium phase. Mathematical models can also be used to provide evidence of such state. The system of partial differential equations proposed in (10) is precisely intended to describe the earliest stages of immune control of tumor growth. Remarkably, while being in the most favorable condition, only taking into account the cytotoxic effector immune cells and no immunosuppressive mechanisms, the model reproduces the formation of an equilibrium phase with maintenance of residual tumor cells rather than their complete elimination. Besides suggesting that elimination may be difficult to reach, this finding also brings out the role of leading parameters that shape the equilibrium features and opens new perspectives to elaborate cancer therapy strategies that reach this state of equilibrium.
To decipher tumor-immune system dynamics leading to equilibrium state, we have developed here computational tools. The total mass of the tumor is a critical criterion of the equilibrium and was used to predict parameters that contribute the most to the establishment of the equilibrium. By means of global sensitivity analysis, we identified one leading parameter, A, and three others, R, γ and a that affect the most the variability of the immune-controlled tumor mass; A, R and γ are related to immune cells, and a to tumor cells. Moreover, the influence of the leading parameters is significantly increased when they are paired. This observation supports the development of combined therapeutic treatments which would be more efficient at reducing tumor growth and toxicity. Because the pairs (A, R), (A, γ), (R, γ), (A, a), (a, R) and (a, γ) are the most influential, we predict that a combination of drugs enhancing antitumor immune responses with drugs diminishing tumor aggressiveness will be the most efficient. This is consistent with the clinical benefit obtained when chemotherapies reducing the tumor cell division rate a are combined with immunotherapies increasing A and R (), The parameter A which governs the efficacy of the immune system to eliminate tumor cells, is the most influential. This finding is consistent with the observation that “hot” tumors infiltrated with immune cells have better prognostic than “cold” tumors () and that the immune cells with the strongest positive impact on patient’s survival are the cytotoxic CD8+ T cells (). It is also in line with the success of ICP which revert immune tolerance triggered by chronic activation and upregulation of exhaustion markers on effector T and NK cells, thus not only increasing the parameter A but also R (). The leading role of the parameter A is also demonstrated by experimental studies and clinical trials, such as adoptive transfer of CAR-T and CAR-NK cells engineered to attack cancer cells, immunomodulating antibody therapies or cancer vaccines which boost the antitumor immune response (, –). Finally, our finding that the parameter γ is highly influential is confirmed by the administration of cytokines that stimulate and increase effector T and NK cell survival which are efficient at controlling tumor growth (). Thus, altogether, these experimental and clinical data validate the numerical method.
Interestingly, besides the dominant role of the parameter A, only two additional parameters related to immune cells R, γ seem to have an influence on the tumor mass at equilibrium. These data predict that to enhance the immune response, it is more efficient to increase the rate of influx and conversion of naive immune cells into effector cells (parameter R) or to increase the lifespan of immune effectors (parameter γ) than to increase chemotaxis as a whole (parameters χ, Aa, K ). The lack of influence of chemotaxis emphasizes that the localization of immune cells within tumors is necessary but not sufficient. Indeed, the leading influence of the parameters A, R, γ stresses the importance of having functional immune cells infiltrating tumors. Overcoming immune suppression is therefore highly relevant in therapeutic strategies.
Targeting Immune-mediated tumor mass dormancy is gaining more and more attention, having been linked to recurrence and metastasis (, ). The persistence of undetectable tumor cells after primary tumor resection at the primary site but also their spreading to metastatic niches are major causes of treatment failure. Thus, developing strategies to maintain an equilibrium between these tumor cells and the immune response is crucial. Interestingly, a recent study demonstrated a role of the NK cell reservoir in blocking the reawakening of dormant tumor cells (). The mechanisms involve IL-15 that drives NK cell proliferation and IFN- γ secreted by NK. Therapies boosting NK cell activity like IL-15 superagonists, or engineered NK cell engagers are therefore promising strategies to sustain NK cell-mediated maintenance of tumor dormancy (, ).
It is appropriate to finally comment on the limitations of this work and provide new avenues for future research. Firstly, the analysis focuses on the asymptotic state, taking full advantage of its mathematical interpretation which makes it easily computable. However, the transient states and the rate at which the equilibrium becomes observable are simply disregarded, while they are certainly essential for assessing the biological relevance of the equilibrium state. Further analysis is therefore needed in order to understand how the parameters of the model influence the trend to equilibrium. Secondly, the modeling approach is facing contradictory requests: on the one hand, the lack of knowledge on the parameters motivates working with a reduced set of equations, at the cost of considering an “averaged” behavior (say for instance between different types of immune cells); on the other hand, it might be important to keep under consideration many relevant and competing effects of cellular interactions. These issues can be addressed with a better access to biological data and through the development of dedicated methods of parameter identification. This is of course even more important when describing the effects of treatments. Thirdly, the present analysis is limited to an idealized situation in which many important effects have been overlooked. In particular, the immune response can also promote the tumor growth. Considering such immune actions leads to a much more complex dynamical behavior and the possible establishment of an escape phase, as shown in (). Finally geometrical aspects and heterogeneity are poorly addressed and restrict the relevance of the description to the earliest stages of the tumor development. More complex models, with a full space structuration, should be elaborated in order to obtain a more accurate description of the tumor microenvironment.
5 Conclusion
In conclusion, clinical trials have been undertaken quite often on assumptions from acquired knowledge on tumor development and immune responses to cancer cells, but without tools to help the decision-making. The numerical methods developed here provide valuable hints for the design and the optimization of antitumor therapies. The approach is in agreement with published experimental findings and clinical evidence. By adapting the range of the parameters to the biological values, one can more precisely adapt the therapeutic strategies to specific types of tumors. We thus conclude that mathematical modelling combined with numerical validation provide valuable information that could contribute to better stratify the patients eligible for treatments and consequently save time and lives. In addition, it could also help to decrease the burden of treatment cost providing hints on optimized therapeutic strategies.
6 Computation of the Eigen-Elements of the Growth-Division Equation
The binary division operator (2) is a very specific case, and for applications it is relevant to deal with more general expressions. Namely, we have
In (), a( z′ ) is the frequency of division of cells having size z′ , and k( z∣z′ ) gives the size-distribution that results from the division of a tumor cell with size z′. What is crucial for modeling purposes is the requirement
which is related to the principle that cell-division does not change the total mass
We refer the reader to () for examples of such cell-division operators and the analysis of the eigenvalue problem (6) under quite general assumptions of the growth rate V, the frequency a and the kernel k. Our numerical method can handle such general coefficients.
It is important to bear in mind the main arguments of the proof of the existence-uniqueness of the eigenpair for the growth-division equation. Namely, for Λ large enough we consider the shifted operator
Then, we check that the operator SΛ which associates to a function f the solution n of TΛn = f fulfills the requirements of the Krein-Rutman theorem (roughly speaking, positivity and compactness), see (). Accordingly, the quantity of interest λ is related to the leading eigenvalue of SΛ. In fact, this reasoning should be applied to a somehow truncated and regularized version of the operator, and the conclusion needs further compactness arguments; nevertheless this is the essence of the proof. In terms of numerical method, this suggests to appeal to the inverse power algorithm, applied to a discretized version of the equation. However, we need to define appropriately the shift parameter Λ. As far as the continuous problem is considered, Λ can be estimated by the parameters of the model (), but it is critical for practical issues to check whether or not this condition is impacted by the discretization procedure. This information will be used to apply the inverse power method to the discretized and shifted version of the problem.
6.1 Analysis of the Discrete Problem
The computational domain for the size variable is the interval [0, R] where R is chosen large enough: due to the division processes, we expect that the support of the solution remains essentially on a bounded interval, and the cut-off should not perturb too much the solution. In what follows, the size step h = zi+1 -zi is assumed to be constant. The discrete unknowns Ni, with i ∈ { 1, … ,I } and h = R/I, are intended to approximate N( zi ) where zi = ih. The integral that defines the gain term of the division operator is approximated by a simple quadrature rule. For the operator (2) the kernel involves Dirac masses which can be approached by peaked Gaussian. We introduce the operator defined by
where Fi =Vi+1/2Ni represents the convective numerical flux on the grid point zi+1/2 =( i+1/2 )h , i ∈ { 1, … ,I } . This definition takes into account that the growth rate is non negative, and applies the upwinding principles. Note that the step size h should be small enough to capture the division of small cells, if any. The following statement provides the a priori estimate which allows us to determine the shift for the discrete problem.
Theorem 1.1. We suppose that
z↦V( z ) is a continuous function which lies in L∞ and it is bounded from below by a positive constant,
remains bounded uniformly with respect to h,
for any i∈{ 1,...,I−1 } , there exists j∈{ 1,...,I−1 } such that a(zj)k( zi ∣zj )>0,
there exists Z0 ∈( 0,∞ ) such that, setting , we have for any z ≥ Z0.
Let
and we suppose that R > Z0 is large enough. Then, is invertible and there exists a pair μ>0, N∈ℝIwith positive components, such that. Moreover.
Note that the sum that defines is actually reduced over the indices such that jh≤z ; this quantity is interpreted as the expected number of cells produced from the division of a cell with size z so that the forth assumption is quite natural.
Proof. Let f∈ℝI. We consider the equation
We denote the solution. We are going to show that is well defined and satisfies the assumptions of the Perron-Frobenius theorem, see e. g. (47, Theorem 1.37 & Corollary 1.39) or (, Chapter 5).
It is convenient to introduce the change of unknown Ui =NiVi+1/2 , ∀i∈{ 1,⋯,I } . The problem recasts as
The solution is interpreted as the fixed point of the mapping
where U is given by U1 = 0 and
We are going to show that Ah is a contraction: ∥Ahξ∥ℓ∞ ≤k∥ξ∥ℓ∞ for some k < 1. Multiplying (20) by sign (Ui), we obtain
We multiply this by the weight , where all factors are ≥1 . We get
Then, summing over i∈{ 2, ... ,m } yields
where actually U1 = 0. It follows that
Therefore, Ah is a contraction provided (19) holds. This estimate is similar to the condition obtained for the continuous problem, see (, Proof of Theorem 2, Appendix B); the discretization does not introduce further constraints.
We are now going to show that is a M-matrix when (19) holds. Let f ∈ ℝI ∖{ 0 } with non negative components. Let U ∈ ℝI satisfy . Let i0 be the index such that Ui0 =min { Ui , i∈{ 2, ..., I } }. We have
Since fi0 ≥0 , we get
which tells us that Ui0 ≥0. Suppose Ui0 =0 for some i0 >1. Coming back to (21), we deduce that Ui0 −1 vanishes too, and so on and so forth, we obtain U1 =⋯Ui0 =0 . Finally, we use the irreductibility assumption iii): we can find j0 > i0 such that and (21) implies , so that Uj0 =0. We deduce that U = 0, which contradicts f≠0 . Therefore the components of U are positive, but U1.
We conclude by applying the Perron-Froebenius theorem to , (). It remains to prove that is positive, with μ the spectral radius of . To this end, we make use of assumption iv). We set Z0 = i0h. We argue by contradiction, supposing that λ=Λ−1/μ<0. We consider the eigenvector with positive components and normalized by the condition . We have
It follows that, for m≥i0 ,
It implies
We arrive at
a contradiction when R is chosen large enough (but how large R should be does not depend on h). Therefore, we conclude that λ > 0.
6.2 Numerical Approximation of (λ, N)
We compute (an approximation of) the eigenpair (λ, N) by using the inverse power method which finds the eigenvalue of with largest modulus:
• We pick Λ verifying (19).
• We compute once for all the LU decomposition of the matrix .
• We choose a threshold 0 < ε ≪ 1.
• We start from a random vector N(0) and we construct the iterations
until the relative error is small enough. Then, given the last iterate N(K), we set
This approach relies on the ability to approximate correctly the eigenpair of the growth-fragmentation operator. In particular, it is important to preserve the algebraic multiplicity. This issue is quite subtle and it is known that the pointwise convergence of the operator is not enough to guarantee the convergence of the eigenelements and the consistency of the invariant subspaces, see () for relevant examples. This question has been thoroughly investigated in (, ) which introduced a suitable notion of stability. It turns out that one needs a uniform convergence of the operators. Namely, here, we should check that as I→∞. In the present framework, a difficulty relies on the fact that the size variable lies in an unbounded domain, which prevents for using usual compactness arguments. For this reason, we introduce a truncated version of the problem, which has also to be suitably regularized. Let us denote by the corresponding operator, where ε represents the regularization parameter. This truncated and regularized perator appeared already in (). Indeed, we know from () that as R→∞ and ε→0, hence, this implies that as R→∞ and ε→0 by continuity of the map Π:TΛ ↦ ( TΛ )−1 . Moreover, is well-defined, continuous and compact, see (32, Appendix. B). The discrete operators converge pointwise to , and the compactness of ensures that the discrete operator converges uniformly to , for 0<R<ε and 0<ε<1 fixed (see () for more details on this fact). Following (), we deduce that the numerical eigenelements (λI, NI) converges to ( λR,ε , NR,ε ), the eigenelements of , while preserving their algebraic multiplicity. Finally the uniform convergence as R→∞ and ∈→0 ensures the convergence of ( λR,ε , NR,ε ), to (λ, N) ().
6.3 Numerical Results
For some specific fragmentation kernels and growth rates, the eigenpair is explicitly known, see (). We can use these formula to check that the algorithm is able to find the expected values and profiles. To this end, we introduce the relative errors
where N(K) and N are both normalized by .
Mitosis fragmentation kernel. We start with the binary division kernel:
The associated division operator is described by (2). We assume that a and V are constant. In this specific case the eigenpair is given by
with N > 0 an appropriate normalizing constant and ( αn )n∈ℕ is the sequence defined by the recursion
In practice we shall use a truncated version of the series that defines N. For the numerical tests, we use the parameters collected in Table 4.
Table 4
| a | V | R | ∈ |
|---|---|---|---|
| 4 | 0.6 | 5 | 10-6 |
Data for the numerical tests: binary division kernel.
With this threshold ε, the approached eigenpair is reached in 43 iterations, independently of the size step. Figure 10 represents the evolution of the error as a function of h in a log-log scale, see also Table 5: N(K) approaches N at order 1. The rate improves when using a quadrature rule with a better accuracy. For this test, the approximation of the eigenvalue is already accurate with a coarse grid; it is simply driven by the threshold ε and does not significantly change with h.
Table 5
| Number of cells | Eλ | EV |
|---|---|---|
| 1000 | 3.73 × 10-5 | 3.83 × 10-2 |
| 2000 | 5.68 × 10-8 | 1.93 × 10-2 |
| 4000 | 6.77 × 10-7 | 9.69 × 10-3 |
| 8000 | 6.84 × 10-7 | 4.85 × 10-3 |
Binary division kernel: errors for several number of grid points.
Figure 10
Uniform fragmentation. The Uniform fragmentation kernel is defined by:
We apply the algorithm for the following two cases:
V(z) = V0 and a(z) = a0z. We have and
We still use the values in Table 4 (especially, a0 = a and V0 = V). The approximated eigenpair is obtained in 84 iterations and, as in the previous test, it does not change with the size step. In this case, both the eigenvalue and the eigenfunction are approached at order 1, see Table 6 and Figure 11.
Table 6
| Number of cells | Eλ | EV |
|---|---|---|
| 1000 | 1.30 × 10-2 | 8.89 × 10-3 |
| 2000 | 6.43 × 10-3 | 4.50 × 10-3 |
| 4000 | 3.23 × 10-3 | 2.24 × 10-3 |
| 8000 | 1.62 × 10-3 | 1.13 × 10-3 |
Uniform fragmentation, ex. 1: errors for several number of grid points.
Figure 11
V(z) = V0z and a(z) = a0zn with n∈ℕ∖{ 0 } . The eigenpair is defined by the following formula:
Note that the growth rate V vanishes and Theorem 1.1 does not apply as such. Nonetheless, the algorithm works well and still captures the eigenpair. We perform the test for n = 1 and n = 2 and the results are recorded in Table 7; Figure 12 and Table 8; Figure 13, respectively.
Table 7
| Number of cells | Eλ | EV |
|---|---|---|
| 1000 | 4.70 × 10-2 | 2 × 10-2 |
| 2000 | 2.43 × 10-2 | 1.06 × 10-2 |
| 4000 | 1.25 × 10-2 | 5.5 × 10-3 |
| 8000 | 6.39 × 10-3 | 2.81 × 10-3 |
Uniform fragmentation, ex. 2, case n = 1: errors for different number of cells.
Figure 12
Table 8
| Number of cells | Eλ | EV |
|---|---|---|
| 1000 | 2.39 × 10-2 | 8.81 × 10-2 |
| 2000 | 1.23 × 10-3 | 4.53 × 10-3 |
| 4000 | 6.41 × 10-3 | 2.35 × 10-3 |
| 8000 | 3.41 × 10-3 | 1.24 × 10-3 |
Uniform fragmentation, ex. 2, case n = 2: errors for different number of cells.
Figure 13
7 Sensitivity Analysis on the Equilibrium Mass
Having an efficient procedure to predict the residual mass of the equilibrium phase also opens perspectives to discuss the influence of the parameters. This can provide useful hints for the design and the optimization of anti-tumor therapies. We address this issue by performing a global sensitivity analysis on the immune-controlled tumor mass. Sensitivity analysis also provides information on the quantification of uncertainty in the model output with respect to the uncertainties in the input parameters. We remind the reader that the equilbrium mass is seen as a function of the parameters in Table 1:
We consider that the input parameters are independent random variables uniformly distributed in an interval [ x1 ,x2 ]⊂( 0,∞ ) :
The pillar of the Sobol sensitivity analysis is the decomposition of f into 2n - 1 summands of increasing dimensions:
Where
and Mi1 ⋯ip =( Mi1 ,⋯Mip ) . The existence and uniqueness of the above decomposition has been proven in (69), given f a square integrable function. Owing to the orthogonality condition (29), the total variance of f reads:
Given (26), V can be decomposed as follows:
where the terms Vi1 ⋯ip , called partial variances read:
Following the description in (), the Sobol’ sensitivity indices are defined as follows:
They verify
Each index Si1 ⋯ip measures how the total variance of f is affected by uncertainties in the set of input parameters i1 ⋯ip. An equivalent definition of the above indices is given by [see ()]:
The total effect of a specific input parameter i is evaluated by the so-called total sensitivity index , the sum of the sensitivity indices which contain i:
where Ci ={ ( i1 ⋯ip ):∃m∈{ 1,...,p }, im =i }. In practice, the sensitivity indices that are needed to discriminate the impact of the parameters are the first, second and total Sobol’ sensitivity indices. The above indices are computed using Monte Carlo simulations combined with a careful sampling of the parameters space in order to reduce the computational load and the number of model evaluations. For this purpose, the following estimators can be derived using two different N samples A and B, see (, ),
Here the notation M−( i1 ⋯ip )l stands for the l-th sample line where we get rid of the points corresponding to the indices i1 ,⋯,ip. The total sensitivity () is given by:
where S-i is the sum of all the sensitivity indices that do not contain the index i. Hence, the total sensitivity index estimator reads:
Where
Funding
This work was supported by the French Government (National Research Agency, ANR) through the “Investments for the Future” programs LABEX SIGNALIFE ANR-11-LABX-0028 and IDEX UCAJedi ANR-15-IDEX-01.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/atsoukevin93/tumorgrowth.
Author contributions
Conception and design: KA, VB, and TG; Development of methodology: KA, VB, and TG; Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): FA, VB, and SK; Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): KA, VB, and TG; Writing, review, and/or revision of the manuscript: KA, FA, VB, and TG; Study supervision: FA, VB, and TG. All authors contributed to the article and approved the submitted version.
Acknowledgments
The authors acknowledge the support of UCAncer, an incentive Université Côte d’Azur network, which has permitted and encouraged this collaboration. We thank the IPMC’s animal house and Imaging/Flow cytometry core facilities for technical support, and Franck Aguila for the artwork of the graphical abstract.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
DunnGPBruceATIkedaHOldLJSchreiberRD. Cancer Immunoediting: From Immunosurveillance to Tumor Escape. Nat Immunol (2002) 3:991–8. doi: 10.1038/ni1102-991
2
DunnGPOldLJSchreiberRD. The Immunobiology Review of Cancer Immunosurveillance and Immunoediting. Immunity (2004) 21:137–48. doi: 10.1016/j.immuni.2004.07.017
3
KoebelCMVermiWSwannJBZerafaNRodigSJOldLJet al. Adaptive Immunity Maintains Occult Cancer in an Equilibrium State. Nature (2007) 450:903–8. doi: 10.1038/nature06309
4
BoonTCouliePGden EyndeBJVvan der BruggenP. Human T-Cell Responses Against Melanoma. Annu Rev Immunol (2006) 24:175–208. doi: 10.1146/annurev.immunol.24.021605.090733
5
DranoffG. Cytokines in Cancer Pathogenesis and Cancer Therapy. Nat Rev Cancer (2004) 4:11–22. doi: 10.1038/nrc1252
6
RabinovichGAGabrilovichDSotomayorEM. Immunosuppressive Strategies That are Mediated by Tumor Cells. Ann Rev Immunol (2007) 25:267–96. doi: 10.1146/annurev.immunol.25.022106.141609
7
SmythMJGodfreyDITrapaniJA. A Fresh Look at Tumor Immunosurveillance and Immunotherapy. Nat Immunol (2001) 2:293–9. doi: 10.1088/0031-9155/57/2/R1
8
WhitesideTL. Immune Suppression in Cancer: Effects on Immune Cells, Mechanisms and Future Therapeutic Intervention. Semin Cancer Biol (2006) 16:3–15. doi: 10.1016/j.semcancer.2005.07.008
9
PhanTGCroucherPI. The Dormant Cancer Cell Life Cycle. Nature (2020) 20:398–411. doi: 10.1038/s41568-020-0263-0
10
AtsouKAnjuèreFBraudVMGoudonT. A Size and Space Structured Model Describing Interactions of Tumor Cells With Immune Cells Reveals Cancer Persistent Equilibrium States in Tumorigenesis. J Theor Biol (2020) 490:110163. doi: 10.1016/j.jtbi.2020.110163
11
KirschnerDEPanettaJC. Modeling Immunotherapy of the Tumor-Immune Interaction. J Math Biol (1998) 37:235–52. doi: 10.1007/s002850050127
12
KonstorumAVellaATAdlerAJLaubenbacherRC. Addressing Current Challenges in Cancer Immunotherapy With Mathematical and Computational Modelling. J R Soc Interface (2017) 14:20170150. doi: 10.1098/rsif.2017.0150
13
LaiXCarsonWEStiffAFriedmanAWesolowskiRDugganM. Modeling Combination Therapy for Breast Cancer With BET and Immune Checkpoint Inhibitors. Proc Nat Acad Sc (2018) 115:5534–9. doi: 10.1073/pnas.1721559115
14
de PillisLGRadunskayaAEWisemanCL. A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth. Cancer Res (2005) 65:7950–8. doi: 10.1158/0008-5472.CAN-05-0564
15
de PillisLRadunskayaA. The Dynamics of an Optimally Controlled Tumor Model: A Case Study. Math Comput Model (2003) 37:1221–44. doi: 10.1016/S0895-7177(03)00133-X
16
KuznetsovVKnottG. Modelling Tumor Regrowth and Immunotherapy. Math Comput Model (2001) 33:1275–87. doi: 10.1016/S0895-7177(00)00314-9
17
LiHWangSXuF. Dynamical Analysis of Tumor-Immune-Help T Cells System. Int J Biomath (2019) 12:1950075. doi: 10.1142/S179352451950075X
18
Robertson-TessiMEl-KarehAGorielyA. A Mathematical Model of Tumor–Immune Interactions. J Theor Biol (2012) 294:56–73. doi: 10.1016/j.jtbi.2011.10.027
19
WilkieKPHahnfeldtP. Modeling the Dichotomy of the Immune Response to Cancer: Cytotoxic Effects and Tumor-Promoting Inflammation. Bull Math Biol (2017) 79:1426–48. doi: 10.1007/s11538-017-0291-4
20
ChaplainMAKuznetsovVAJamesZHDavidsonFStepanovaLA. Spatio-Temporal Dynamics of the Immune System Response to Cancer. In: HornMAGieriSWebbGF, editors. Mathematical Models in Medical and Health Science. Vanderbilt, TN, USA: Vanderbilt University Press (1998). International Conference on Mathematical Models in Medical and Health Sciences.
21
EftimieR. Investigation Into the Role of Macrophages Heterogeneity on Solid Tumour Aggregations. Math Biosc (2020) 322:108325. doi: 10.1016/j.mbs.2020.108325
22
KuznetsovVAMakalkinIATaylorMAPerelsonAS. Nonlinear Dynamics of Immunogenic Tumors: Parameter Estimation and Global Bifurcation Analysis. Bull Math Biol (1994) 56:295–321. doi: 10.1007/BF02460644
23
MatzavinosAChaplainMAJKuznetsovVA. Mathematical Modelling of the Spatio-Temporal Response of Cytotoxic T-Lymphocytes to a Solid Tumour. Math Med Biol (2004) 21:1–34. doi: 10.1093/imammb/21.1.1
24
de PillisLGMalletDGRadunskayaAE. Spatial Tumor-Immune Modeling. Comput Math Methods Med (2006) 7:159–76. doi: 10.1080/10273660600968978
25
EftimieRBramsonJLEarnDJD. Interactions Between the Immune System and Cancer: A Brief Review of non-Spatial Mathematical Models. Bull Math Biol (2011) 73:2–32. doi: 10.1007/s11538-010-9526-3
26
EftimieRGillardJJCantrellDA. Mathematical Models for Immunology: Current State of the Art and Future Research Directions. Bull Math Biol (2017) 78:2091–134. doi: 10.1007/s11538-016-0214-9
27
EladdadiAKimPMalletD. Mathematical Models of Tumor-Immune System Dynamics, Springer Proceedings in Math. & Statistics Vol. 107. SwitzerlandSpringer (2014).
28
MahlbacherGEReihmerKCFrieboesaHB. Mathematical Modeling of Tumor-Immune Cell Interactions. J Theor Biol (2019) 469:47–60. doi: 10.1016/j.jtbi.2019.03.002
29
RooseTChapmanSJMainiPK. Mathematical Models of Avascular Tumor Growth. SIAM Rev (2007) 49:179–208. doi: 10.1137/S0036144504446291
30
BellouquidADelitalaM. Mathematical Methods and Tools of Kinetic Theory Towards Modelling Complex Biological Systems. Math Mod Meth Appl Sci (2005) 15:1639–66. doi: 10.1142/S0218202505000923
31
Bekkal BrikciFClairambaultJRibbaBPerthameB. An Age-and-Cyclin-Structured Cell Population Model for Healthy and Tumoral Tissues. J Math Biol (2008) 57:91–110. doi: 10.1007/s00285-007-0147-x
32
Doumic-JauffretMGabrielP. Eigenelements of a General Aggregation-Fragmentation Model. Math Models Methods Appl Sci (2010) 20:757–83. doi: 10.1142/S021820251000443X
33
MichelPMischlerSPerthameB. General Relative Entropy Inequality: An Illustration on Growth Models. J Math Pures Appl (2005) 84:1235–60. doi: 10.1016/j.matpur.2005.04.001
34
MichelP. Existence of a Solution to the Cell Division Eigenproblem. Models Math Meth App Sci (2006) 16:1125–53. doi: 10.1142/S0218202506001480
35
MichelP. Principe D’entropie Généralisée Et Dynamique De Populations Structurées. Ph.D. Thesis, Paris Dauphine. (2005).
36
PerthameB. Transport Equations in Biology. Front Math (Birkhauser) (2007). doi: 10.1007/978-3-7643-7842-4
37
PreziosiL. Modeling Cancer Growth. Chapman Hall: CRC-Press (2003).
38
AlzahraniEOKuangY. Nutrient Limitations as an Explanation of Gompertzian Tumor Growth. Disc Cont Dyn Syst B (2016) 21:357–72. doi: 10.3934/dcdsb.2016.21.357
39
GyllenbergMWebbGF. A Nonlinear Structured Population Model of Tumor Growth With Quiescence. J.Math Biol (1990) 28:671–94. doi: 10.1007/BF00160231
40
IwataKKawasakiKShigesadaN. A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors. J Theor Biol (2000) 203:177–86. doi: 10.1006/jtbi.2000.1075
41
BenzekrySLamontCBeheshtiATraczAEbosJHlatkyLet al. Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth. PloS Comput Biol (2014) 10:e1003800. doi: 10.1371/journal.pcbi.1003800
42
AtsouKAnjuèreFBraudVMGoudonT. A Size and Space Structured Model of Tumor Growth Describes a Key Role for Protumor Immune Cells in Breaking Equilibrium States in Tumorigenesis. PloS One (2021), 0259291. doi: 10.1371/journal.pone.0259291
43
BaccelliFMcDonaldDReynierJ. A Mean Field Model for Multiple TCP Connections Through a Buffer Implementing RED. Perform Eval (2002) 49:77–97. doi: 10.1016/S0166-5316(02)00136-0
44
PerthameBRyzhikL. Exponential Decay for the Fragmentation or Cell-Division Equation. J Differ Eq (2005) 210:155–77. doi: 10.1016/j.jde.2004.10.018
45
FagetJGroeneveldSBoivinGSankarMZanggerNGarciaMet al. Neutrophils and Snail Orchestrate the Establishment of a Pro-Tumor Microenvironment in Lung Cancer. Cell Rep (2017) 21:3190–204. doi: 10.1016/j.celrep.2017.11.052
46
LiCYShanSHuangQBraunRDLanzenJHuKet al. Initial Stages of Tumor Cell-Induced Angiogenesis: Evaluation via Skin Window Chambers in Rodent Models. J Nat Cancer Inst (2000) 92:143–7. doi: 10.1093/jnci/92.2.143
47
GoudonT. Mathematics for Modeling and Scientific Computing. Math Stat (Wiley-ISTE) (2016) 27(3):665–679. doi: 10.1002/9781119371120
48
ChatelinF. Convergence of Approximation Methods to Compute Eigenelements of Linear Operations. SIAM J Numer Anal (1973) 10:939–48. doi: 10.1137/0710080
49
ChatelinF. The Spectral Approximation of Linear Operators With Applications to the Computation of Eigenelements of Differential and Integral Operators. SIAM Rev (1981) 23:495–522. doi: 10.1137/1023099
50
FarrellBEDanieleRPLauffenburgerDA. Quantitative Relationships Between Single-Cell and Cell-Population Model Parameters for Chemosensory Migration Responses of Alveolar Macrophages to C5a. Cell Motil (1990) 16:279–93. doi: 10.1002/cm.970160407
51
FriedmanAHaoW. The Role of Exosomes in Pancreatic Cancer Microenvironment. Bull Math Biol (2018) 80:1111–33. doi: 10.1007/s11538-017-0254-9
52
YatesACallardR. Cell Death and the Maintenance of Immunological Memory. Disc Cont Dyn Syst.-B (2001) 1:43–59. doi: 10.3934/dcdsb.2001.1.43
53
BeckRJSlagterMBeltmanJB. Contact-Dependent Killing by CytotoxicTLymphocytes is Insufficient for EL4 Tumor Regression In Vivo. Cancer Res (2019) 79:3406–16. doi: 10.1158/0008-5472.CAN-18-3147
54
CazauxMGrandjeanCLLemaîtreFGarciaZBeckRJMiloIet al. Single-Cell Imaging of CAR T-Cell Activity In Vivo Reveals Extensive Functional and Anatomical Heterogeneity. J Exp Med (2019) 216:1038–49. doi: 10.1084/jem.20182375
55
HansonHLDonermeyerDLIkedaHWhiteJMShankaranVOldLJet al. Eradication of Established Tumors by CD8+ T-Cell Adoptive Immonutherapy. Cell Press (2000) 13:265–76. doi: 10.1016/s1074-7613(00)00026-1
56
NolzJCHillAB. Strength in Numbers: Visualizing CTL-Mediated Killing In Vivo. Immunity (2016) 44:207–8. doi: 10.1016/j.immuni.2016.01.026
57
KwokCSColeSELiaoSK. Uptake Kinetics of Monoclonal Antibodies by Human Malignant Melanoma Multicell Spheroids. Cancer Res (1988) 48:1856–63.
58
CairnsCMGordonJRLiFBaca-EstradaMEMoyanaTXiangJ. Lymphotactin Expression by Engineered Myeloma Cells Drives Tumor Regression: Mediation by CD4+ and CD8+ T-cells and neutrophils expressing XCR1 receptor. J Immun (2001) 167:57–65. doi: 10.4049/jimmunol.167.1.57
59
KhouSPopaALuciCBihlFMeghraoui-KheddarABourdelyPet al. Tumor-Associated Neutrophils Dampen Adaptive Immunity and Promote Cutaneous Squamous Cell Carcinoma Development. Cancers (2020) 10-12:E1860. doi: 10.3390/cancers12071860
60
SamsonALavielleMMentréF. Extension of the SAEM Algorithm to Left-Censored Data in Nonlinear Mixed-Effects Model: Application to HIV Dynamics Model. Comput Stat Data Anal (2006) 51:1562–74. doi: 10.1016/j.csda.2006.05.007
61
[Dataset] Monolix Version 2020R1. Antony, france: Lixoft sas (2020).
62
WeiseKPoßnerLMüllerEGastRKnöscheTR. Pygpc: A Sensitivity and Uncertainty Analysis Toolbox for Python. SoftwareX (2020) 11:100450. doi: 10.1016/j.softx.2020.100450
63
HermanJUsherW. SALib: An Open-Source Python Library for Sensitivity Analysis. J Open Source Softw (2017) 2:joss00097. doi: 10.21105/joss.00097
64
ShashniBAriyasuSTakedaRSuzukiTShiinaSAkimotoKet al. Size-Based Differentiation of Cancer and Normal Cells by a Particle Size Analyzer Assisted by a Cell-Recognition Pc Software. Biol Pharm Bull (2018) 41:487–503. doi: 10.1248/bpb.b17-00776
65
MooreGESandbergAAWatneAL. The Comparative Size and Structure of Tumor Cells and Clumps in the Blood, Bone Marrow, and Tumor Imprints. Cancer (1960) 13:111–7. doi: 10.1002/1097-0142(196001/02)13:1<111::aid-cncr2820130121>3.0.co;2-y
66
ErdiYE. Limits of Tumor Detectability in Nuclear Medicine and PET. Mol Imaging Radionuclid Ther (2012) 21:23–8. doi: 10.4274/Mirt.138
67
HarenbergDMarelliSSudretBWinschelV. Uncertainty Quantification and Global Sensitivity Analysis for Economic Models. Quanti Econom (2019) 10:1–41. doi: 10.3982/qe866
68
SaltelliA. Sensitivity Analysis of Model Output. An Investigation of New Techniques. Comput Stat Data Anal (1993) 15:211–38. doi: 10.1016/0167-9473(93)90193-W
69
Sobol’I. Sensitivity Estimates for Nonlinear Mathematical Models. Math Modell Comput (1993) 1:407–14.
70
CrestauxTLe MaîtreOMartinezJM. Polynomial Chaos Expansion for Sensitivity Analysis. Reliabil Eng Syst Saf (2009) 94:1161–72. doi: 10.1016/j.ress.2008.10.008
71
SudretB. Global Sensitivity Analysis Using Polynomial Chaos Expansions. Reliabil Eng Syst Saf (2008) 93:964–79. doi: 10.1016/j.ress.2007.04.002
72
ErnstOGMuglerAStarkloffHJUllmannE. On the Convergence of Generalized Polynomial Chaos Expansions. ESAIM: Math Model Numeric Anal (2012) 46:317–39. doi: 10.1051/m2an/2011045
73
Le GratietLMarelliSSudretB. Metamodel-Based Sensitivity Analysis: Polynomial Chaos Expansions and Gaussian Processes. In: GhanemRHigdonDOwhadiH, editors. Handbook of Uncertainty Quantification. Cham: Springer International Publishing (2017). p. 1289–325. doi: 10.1007/978-3-319-12385-1_38
74
SaltelliA. Making Best Use of Model Evaluations to Compute Sensitivity Indices. Comput Phys Commun (2002) 145:280–297. doi: 10.1016/S0010-4655(02)00280-1
75
BaillyCThuruXQuesnelB. Combined Cytotoxic Chemotherapy and Immunotherapy of Cancer: Modern Times. NAR Cancer (2020) 2:1–20. doi: 10.1093/narcan/zcaa002
76
GalonJBruniD. Approaches to Treat Immune Hot, Altered and Cold Tumours With Combination Immunotherapies. Nat Rev Drug Discovery (2019) 18:197–218. doi: 10.1038/s41573-018-0007-y
77
GalonJCostesASanchez-CaboFKirilovskyAMlecnikBLagorce-PagèsCet al. Type, Density, and Location of Immune Cells Within Human Colorectal Tumors Predict Clinical Outcome. Science (2006) 313:1960–4. doi: 10.1126/science.1129139
78
SharmaPAllisonJ. The Future of Immune Checkpoint Therapy. Science (2015) 348:56–61. doi: 10.1126/science.aaa8172
79
ChampiatSTselikasLFarhaneSRaoultTTexierMLanoyEet al. Intratumoral Immunotherapy: From Trial Design to Clinical Practice. Clin Cancer Res (2020) 27(3):665–79. doi: 10.1158/1078-0432.CCR-20-0473
80
KaufmanHLAtkinsMBSubediPWuJChambersJMattinglyTJet al. The Promise of Immuno-Oncology: Implications for Defining the Value of Cancer Treatment. J ImmunoTher Cancer (2019) 7:129–40. doi: 10.1186/s40425-019-0594-0
81
ShekarianTValsesia-WittmannSCauxCMarabelleA. Paradigm Shift in Oncology: Targeting the Immune System Rather Than Cancer Cells. Mutagenesis (2015) 30:205–11. doi: 10.1093/mutage/geu073
82
WangHWangSHuangMLiangXTangYTangY. Targeting Immune-Mediated Dormancy: A Promising Treatment of Cancer. Front Oncol (2019) 9:498:1–9. doi: 10.3389/fonc.2019.00498
83
CorreiaALGuimaraesJCder MaurPASilvaDDTrefnyMPOkamotoRet al. Hepatic Stellate Cells Suppress NK Cell-Sustained Breast Cancer Dormancy. Nature (2021) 594:566–92. doi: 10.1038/s41586-021-03614-z
84
LopesNVivierE. Natural Killer Cells Lull Tumours Into Dormancy. Nature (2021) 594:501–2. doi: 10.1038/d41586-021-01381-5
85
KreinMRutmanM. Linear Operator Leaving Invariant a Cone in a Banach Space. Amer Math Soc Transl (1962) 10:199–325.
86
SerreD. Matrices: Theory and Applications. In: Graduate Texts in Math, vol. 216. New-York, Dordrecht, Heidelberg, London: Springer (2002).
87
HommaTSaltelliA. Importance Measures in Global Sensitivity Analysis of Nonlinear Models. Reliabil Eng Syst Saf (1996) 52:1–17. doi: 10.1016/0951-8320(96)00002-6
Summary
Keywords
cancer, mathematical oncology, equilibrium phase, immunotherapy, drug response
Citation
Atsou K, Khou S, Anjuère F, Braud VM and Goudon T (2022) Analysis of the Equilibrium Phase in Immune-Controlled Tumors Provides Hints for Designing Better Strategies for Cancer Treatment. Front. Oncol. 12:878827. doi: 10.3389/fonc.2022.878827
Received
18 February 2022
Accepted
02 May 2022
Published
27 June 2022
Volume
12 - 2022
Edited by
George Bebis, University of Nevada, Reno, United States
Reviewed by
Maria Ramirez, University of Texas MD Anderson Cancer Center, United States; Caterina AM La Porta, University of Milan, Italy
Updates
Copyright
© 2022 Atsou, Khou, Anjuère, Braud and Goudon.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Véronique M. Braud, braud@ipmc.cnrs.fr; Thierry Goudon, thierry.goudon@inria.fr
‡These authors have contributed equally to this work and share last authorship
†Present address: Sokchea Khou, Department of Cell, Developmental & Cancer Biology, Oregon, Health & Science University, Portland, OR, United States
This article was submitted to Cancer Immunity and Immunotherapy, a section of the journal Frontiers in Oncology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.