Analysis of the Equilibrium Phase in Immune-Controlled Tumors Provides Hints for Designing Better Strategies for Cancer Treatment

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.


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 (1)(2)(3): 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 (4)(5)(6)(7)(8). 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 (9). The equilibrium state is however difficult to apprehend experimentally because the tumor mass at equilibrium is below detectable limits (3). Mathematical modeling of the tumorimmune 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 (10)(11)(12)(13).
We previously (10) introduced a specific multiscale mathematical model based on partial differential equations (PDE), intended to describe the earliest stages of tumorimmune 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 (11,(14)(15)(16)(17)(18)(19) Extension to the PDE framework has permitted to bring out the role of space organisation (20)(21)(22)(23). 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 (24)(25)(26)(27)(28)(29). The original model developed in (10) 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 (10) 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 equilibriumnamely, the size distribution of the tumor cells and the residual tumor massfor 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.

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 (10), summarized by Figure 1, led to couple an evolution equation for the size-distribution of the tumor cells, and a convectiondiffusion 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 W (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 ∈ W .
The unknowns are • The size density of tumor cells ( t, z )↦n( t, z ) so that the integral Z b a z n(t, z) dz 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 ) ↦f( t, x ).
The specific biological assumptions made to construct the model are fully described in (10). Figure 2 offers an overview of the interaction mechanisms embodied in the equations and of the role of the parameters of the model. 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 zn t, z ð Þ dz : 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) FIGURE 1 | Schematic view of the geometry of the mathematical model. The tumor cells are located at the center of the domain where they are subjected to growth and division mechanisms. Immune cells are activated from baths of resting cells; their motion is driven by diffusion combined to a convection field, due to chemotactic mechanisms and directed towards the tumor. growth rate z↦V( z )≥0 and the cell division operator Q(n). We refer the reader to (30)(31)(32)(33)(34)(35)(36)(37) 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 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 (32) (see Appendix A). The specific case where the division rate a in (2) is a positive constant makes the model simpler, and is often used. It is however likely relevant to incorporate more complex behaviors through the sizedependence; 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 (38)(39)(40); a detailed study of growth laws can be found in (41). 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 cellsat 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 equationthe death rate g>0, the chemotactic strength c>0 , and the diffusion coefficient Dcorrespond 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 m 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 f, 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↦s( 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 ∂W. Finally, the activated immune cells are able to destroy tumor cells, as described by the death term in (1a) where the positive parameters A,A s and q,q s 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 (10) 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 (42) on these issues. Remarkably, the model (1a)-(1e) is able to reproduce equilibrium phases where the tumor growth is controlled by the immune response.

Development of Numerical Methods Predicting Parameters of the Equilibrium in Immune-Controlled Tumors
According to (2,3,9), 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 (10) 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.

Equilibrium States
The definition of the equilibrium relies on the following arguments. When disregarding the immune response, the celldivision equation admits a positive eigenstate, which drives the large time behavior of the solution. To be more specific, there exists l>0 and a non negative function z ≥ 0 ↦ N(z) satisfying The existence-uniqueness of the eigenpair (l, N) can be found in (32,34). Furthermore, when the tumor does not interact with the immune system, the large time behavior is precisely driven by the eigenpair: the solution of (5) behaves like where m 0 >0 is a constant determined by the initial condition, see (33,34). Consequently, in the immune-free case, the tumor population grows exponentially fast, with a rate l>0 , and, as time becomes large, its size repartition obeys a certain profile N. In the specific case where V is constant and Q is the binary division operator (2), with a constant division rate a, we simply have l=a and the profile N is explicitly known (43,44). 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 (l, N). Coming back to the coupled model (1a)-(1e), we infer that the equilibrium phase corresponds to the situation where the death ratethe integral of the immune cells concentration with weight d, denoted as m c 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 m 0 N(z). The proportionnality factor is related to the total mass by the relation m 1 = m 0 • 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 Z This can be interpreted as an implicit definition of the total mass m 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 m 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 )∈L 2 ( W ) be a non negative function. If l>0 is small enough, there exists a unique m 1 ( l )>0 such that the solution C m 1 ( l ) 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 m c (t) ∼ l as time becomes large, and the size repartition of tumor cells tends to the eigenfunction N. The total mass m 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 (45). 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 (45,46). 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 (l, N) 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.

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: l 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 [(47), Section 1.2.5]. We refer the reader to (48,49) 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 (32), 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 (l, N) is explicitly known, see (32). We used these formula to validate the ability of the algorithm to find the expected values and profiles.

Computation of the Equilibrium Mass
Having at hand the eigenvalue l, we go back to the convectiondiffusion equation (7) and the constraint (8) that determine implicitly the total mass m 1 of the residual tumor. For a given value of m 1 , we numerically solve (7) by using a finite volume scheme, see (10, 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 < m a < m b ; the mass we are searching for is expected to belong to the interval (m a , m b ). • Set m 1 = m a +m b 2 and compute the associated solution C m 1 of (7) (the subscript emphasizes the dependence with respect to m 1 ). Evaluate the discrete version of the quantity I=∫ dC m 1 dx−l • If I < 0, then replace m a by m 1 , otherwise replace m b by m 1 .
• We stop the algorithm when the relative error m b −m a m a < e 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.

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 (59) where the development of chemically-induced cutaneous squamous cell carcinoma (cSCC) is investigated.
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 (32)(33)(34)44). Let t ↦ m 0 (t) = 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 d dt m 0 = am 0 , d dt m 1 = Vm 0 : 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 ij is a typical observation of the mouse i for a given measurement type k ∈ { 0,1,2} (with (0, 1, 2) referring to ( m 0 ,m 1 , c ) respectively) at time t k ij for i ∈ { 1, … ,N } and i ∈ { 1, … ,N } . We suppose that the statistics of the measurements obeys, for where f (k) (t k ij ; q k i ) is the evaluation of the model at time t k ij , q k i ∈ R p is the vector of the parameters describing the individual i and e (k) ij 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 h k i which explain the interindividual 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( w 0 ,w 1 ,w 2 ) where the w 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.

Estimation of the Model Parameters
According to the experimental procedure in (59), 5× 10 5 mSCC38 were injected to each mouse at time t 0 = 0. Therefore we fixed the initial number of tumor cells to m 0 ( 0 )=5× 10 5 cells.
Assuming that each tumor cell is spherically shaped with a radius 15 m m, we set m 1 ( 0 )=7.1mm 3 . The initial concentration of immune cells is fixed to c 0 = 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 (60). Let I k ij be the finite or infinite censoring interval for mouse i, measurement k and time t k ij and 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 (61). Furthermore, the individual parameter estimatorsq k i are computed in MONOLIX (61) by means of the Empirical Bayes Estimate (EBE) of q k i which corresponds to the mode of the conditional distribution p(q k i | y k i ;â ) (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 )=r aR w a w R , where r aR represents the correlation coefficient between a and R, enhances the goodness of fit. The estimated value of r 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 h (k) i (the p-values for h a ,h V and h 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.

Mice
FVB/N wild-type (WT) mice (Charles River Laboratories, St Germain Nuelles, France) were bred and housed in specificpathogen-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 Cote d'Azur, France).

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 (61), the library Pygpc for the generalized Polynomial Chaos approximation (62) and the library Salib for the sensitivity analysis (63).

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 W, 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 m 0 (0) = 1 cell n , ue m 1 (0) = 14137.2 mm 3 , 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 tends to the eigenvalue of the cell-division equation, the total mass m 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 mm, which is consistent with observational data reporting the mean size distribution of cancer cells (64,65). For the simplest model of growth-division with a and V constant, we know an expression of the eigenstate (l, N); 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 m f 1 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 m pd 1 be the mass predicted by the powerdichotomy 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.
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 (42).

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 It is important to note that the predicted diameter of the tumor at equilibriumsee Figure 5 is significantly below modern clinical PET scanners resolution limit, which could detect tumors with a diameter larger than 7mm (66). 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 g reduces the efficacy of the immune response and increases the equilibrium tumor size see Figure 6D.
Moreover, as mentioned above, not only the parameters determine the equilibrium mass, but they also impact how the equilibrium establishes.  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.

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 c, the death rate g of the immune cells, and the strength of the chemical signal induced by each tumor cell A s • 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 (67), 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,c,D,A s ,g,K ) and by lognormal distributions for the parameters a and R. In what follows, the total mass at equilibrium, m 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 m 1 of the algorithm is influenced by some subsets i 1 ⋯i p of the input parameters i 1 ⋯i p ( 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 S (i) T , the sum of the sensitivity indices which contain the parameter i.
The horizontal dotted line represents the predicted asymptotic value for m c . The solid line gives the envelope of the oscillations, indicating a polynomial damping rate. The equilibrium needs more time to establish as the strength of the immune system decreases.
can be found in Appendix B). The computation of these indices is usually based on a Monte Carlo (MC) method [see (68,69)] which requires a large number of evaluations of the model due to its slow convergence rate (O(1= ffiffiffiffi N p ) 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 (70)]. 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 q( w )=( a,A,R,c,D,A s ,g,K ), where w represents an element of the set of possible outcomes. Further details on the method can be found in (71). For uniform distributions, the most suitable orthonomal polynomial basis is the Legendre polynomials. The analysis of the distribution of m 1 after a suitable sampling of the parameters space indicates that m 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 (72)]. 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 m 1 . Typically, ln(m 1 ) is decomposed as follows: where ϵ corresponds to the approximation error, J k,p = fa ∈ N k : o k i=1 a i ≤ pg and p represents the highest degree of the expansion. Hence, the dimension of the polynomial basis is given by (k+p) ! k ! p ! . 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 (73). The resulting LOO error indicates 0.4% prediction error. The Sobol's sensitivity indices are then computed from the exponential of the surrogate model (16) by using Monte Carlo simulations combined with a careful space-filling sampling of the parameters space [see (68,74)]. For the computations, a sample with N=1.8× 10 6 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.
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.

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 g 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 g. 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 c, like the strength of the chemical signal induced by each tumor cell A s , 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,g), (R, g), (a,A), (a, R) and (a, g). 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 g 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.

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, g. 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 (75).

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

Atsou et al.
Equilibrium Phase in Immune-Controlled Tumors 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 g; and the tumor cell division rate a with the death rate g. 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

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 tumorrelated 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 g) 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, g and a that affect the most the variability of the immune-controlled tumor mass; A, R and g 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, g), (R, g), (A, a), (a, R) and (a, g) 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 (75), 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 (76) and that the immune cells with the strongest positive impact on patient's survival are the cytotoxic CD8 + T cells (77). 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 (78). 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 (75,(79)(80)(81). Finally, our finding that the parameter g 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 (81). 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, g 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 g) than to increase chemotaxis as a whole (parameters c, A a , 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, g 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 (9,82). 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 (83). The mechanisms involve IL-15 that drives NK cell proliferation and IFN-g 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 (83,84).
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 (42). 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.

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.

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 (21), 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 Z ∞ 0 zQ n ð Þ dz = 0: We refer the reader to (32) 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 (l, N) for the growth-division equation. Namely, for L large enough we consider the shifted operator Then, we check that the operator S L which associates to a function f the solution n of T L n = f fulfills the requirements of the Krein-Rutman theorem (roughly speaking, positivity and compactness), see (85). Accordingly, the quantity of interest l is related to the leading eigenvalue of S L . 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 L. As far as the continuous problem is considered, L can be estimated by the parameters of the model (32), 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.

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 = z i+1 -z i is assumed to be constant. The discrete unknowns N i , with i ∈ { 1, … ,I } and h = R/I, are intended to approximate N( z i ) where z i = 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 T h L : where F i =V i+1/2 N i represents the convective numerical flux on the grid point z i+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. ii. h o I j=1 a(z j )k(z i | z j ) remains bounded uniformly with respect to h, iii. for any i∈{ 1,...,I−1 } , there exists j∈{ 1,...,I−1 } such that a(z j )k( z i |z j )>0, iv. there exists Z 0 ∈( 0,∞ ) such that, setting and we suppose that R > Z 0 is large enough. Then, T h L is invertible and there exists a pair m>0, N∈R I with positive components, such thatKer((T h L ) −1 − m) = SpanfNg. Moreover l = L − 1 m > 0. Note that the sum that defines N (z) 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∈R I . We consider the equation It is convenient to introduce the change of unknown U i =N i V i+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 U 1 = 0 and We are going to show that A h is a contraction: ∥A h x∥ ℓ ∞ ≤k∥x∥ ℓ ∞ for some k < 1. Multiplying (20) by sign (U i ), we obtain We multiply this by the weight , where all factors are ≥1 . We get Then, summing over i∈{ 2, ... ,m } yields where actually U 1 = 0. It follows that Therefore, A h is a contraction provided (19) holds. This estimate is similar to the condition obtained for the continuous problem, see (32, Proof of Theorem 2, Appendix B); the discretization does not introduce further constraints.
We are now going to show that T h L is a M-matrix when (19) holds. Let f ∈ R I ∖{ 0 } with non negative components. Let U ∈ R I satisfy (T h L U) i = h f i V i+1=2 . Let i 0 be the index such that U i 0 =min { U i , i∈{ 2, ..., I } }. We have Since f i 0 ≥0 , we get Coming back to (21), we deduce that U i 0 −1 vanishes too, and so on and so forth, we obtain U 1 =⋯U i 0 =0 . Finally, we use the irreductibility assumption iii): we can find j 0 > i 0 such that (21) implies . We deduce that U = 0, which contradicts f≠0 . Therefore the components of U are positive, but U 1 .
We conclude by applying the Perron-Froebenius theorem to (T h L ) −1 , (86). It remains to prove that l = L − 1 m is positive, with m the spectral radius of (T h L ) −1 . To this end, we make use of assumption iv). We set Z 0 = i 0 h. We argue by contradiction, supposing that l=L−1/m<0. We consider the eigenvector with positive components and normalized by the condition h o We have It follows that, for m≥i 0 , a contradiction when R is chosen large enough (but how large R should be does not depend on h). Therefore, we conclude that l > 0.

Numerical Approximation of (l, N)
We compute (an approximation of) the eigenpair (l, N) by using the inverse power method which finds the eigenvalue of (T h L ) −1 with largest modulus: • We pick L verifying (19).
• We compute once for all the LU decomposition of the matrix T h L . • We choose a threshold 0 < e ≪ 1.
• We start from a random vector N (0) and we construct the iterations , and l = L − 1=m : 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 (48) for relevant examples. This question has been thoroughly investigated in (48,49) 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 ∥ (T I L ) −1 − (T L ) −1 ∥ ! 0 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 T R,e L the corresponding operator, where e represents the regularization parameter. This truncated and regularized perator appeared already in (32). Indeed, we know from (32) that ∥ T R,e L − T L ∥ ! 0 as R!∞ and e!0, hence, this implies that ∥ (T R,e L ) −1 − (T L ) −1 ∥ ! 0 as R!∞ and e!0 by continuity of the map P:T L ↦ ( T L ) −1 . Moreover, (T R,e L ) −1 is well-defined, continuous and compact, see (32, Appendix. B). The discrete operators (T I L ) −1 converge pointwise to (T R,e L ) −1 , and the compactness of (T R,e L ) −1 ensures that the discrete operator converges uniformly to (T R,e L ) −1 , for 0<R<e and 0<e<1 fixed (see (49) for more details on this fact). Following (49), we deduce that the numerical eigenelements (l I , N I ) converges to ( l R,e , N R,e ), the eigenelements of (T R,e L ) −1 , while preserving their algebraic multiplicity. Finally the uniform convergence ∥ (T R,e L ) −1 − (T L ) −1 ∥ ! 0 as R!∞ and ∈!0 ensures the convergence of ( l R,e , N R,e ), to (l, N) (32).

Numerical Results
For some specific fragmentation kernels and growth rates, the eigenpair (l, N) is explicitly known, see (32). 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 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 ( a n ) n∈N is the sequence defined by the recursion a 0 = 1, a n = 2 2 n −1 a n−1: 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.
With this threshold e, the approached eigenpair is reached in 43 iterations, independently of the size step. Figure 10 represents the evolution of the error E h V 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 e and E h L does not significantly change with h.
Uniform fragmentation. The Uniform fragmentation kernel is defined by: k(z | z 0 ) = 1 z 0 1 0≤z≤z 0 : We apply the algorithm for the following two cases: We still use the values in Table 4 (especially, a 0 = a and V 0 = 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. V(z) = V 0 z and a(z) = a 0 z n with n∈N∖{ 0 } . The eigenpair is defined by the following formula: exp − a 0 nV 0 z n 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 FIGURE 10 | Binary division kernel: convergence rates of (l (K) , N (K) ) with respect to h.  and the results are recorded in Table 7; Figure 12 and Table 8; Figure 13, respectively.

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 [ x 1 ,x 2 ]⊂( 0,∞ ) : The pillar of the Sobol sensitivity analysis is the decomposition of f into 2 n -1 summands of increasing dimensions: Where 1 x 2 −x 1 Z x 1 ,x 2 ½ f i 1 ⋯ i p M i 1 ⋯ i p dM i k = 0 for k ∈ 1, :::, p f g, (27) f 0 = 1 and M i 1 ⋯i p =( M i 1 ,⋯M i p ) . 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 V i 1 ⋯i p , called partial variances read: Following the description in (69), the Sobol' sensitivity indices are defined as follows: Each index S i 1 ⋯i p measures how the total variance of f is affected by uncertainties in the set of input parameters i 1 ⋯i p . An equivalent definition of the above indices is given by [see (68)]: The total effect of a specific input parameter i is evaluated by the so-called total sensitivity index S (i) T , the sum of the sensitivity indices which contain i: where C i ={ ( i 1 ⋯i p ):∃m∈{ 1,...,p }, i m =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 (68,74), Here the notation M −( i 1 ⋯i p )l stands for the l-th sample line where we get rid of the points corresponding to the indices i 1 ,⋯, i p . The total sensitivity (87) 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: Wherê

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.

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.