Bayesian inference of a non-local proliferation model

From a systems biology perspective, the majority of cancer models, although interesting and providing a qualitative explanation of some problems, have a major disadvantage in that they usually miss a genuine connection with experimental data. Having this in mind, in this paper, we aim at contributing to the improvement of many cancer models which contain a proliferation term. To this end, we propose a new non-local model of cell proliferation. We select data that are suitable to perform Bayesian inference for unknown parameters and we provide a discussion on the range of applicability of the model. Furthermore, we provide proof of the stability of posterior distributions in total variation norm which exploits the theory of spaces of measures equipped with the weighted flat norm. In a companion paper, we provide detailed proof of the well-posedness of the problem and we investigate the convergence of the escalator boxcar train (EBT) algorithm applied to solve the equation.


Introduction
Mathematical models of complex biological phenomena that are developed nowadays are based on the knowledge of biophysical processes. Thanks to this, we correctly obtain the general structure of equations, but unfortunately, we miss the information on parameter values. This is obviously mainly due to the difficulty of performing specific experimental measurements. It is important to underline that sometimes there are no established experimental scenarios with which one can measure the desired parameter values to say nothing about the errors related to measurements when possible.
The lack of reliable values of the parameters implies a significant limitation of the applicability of those models. To overcome this difficulty the approach taken so far was, to search for the model parameters in the literature. However, these were usually deficient or obtained within particular experimental regimes, very often not corresponding to the scenario under consideration. Therefore it is not clear at all that their interpretation is like the parameters in the mathematical model. Many, if not most, cancer models contain some type of logistic function to describe cancer cell proliferation, see [2,4,33,34,44] and references therein. Mainly due to its conceptual ease this short of reliable values of the parameters. Within this paper, we aim at improving potential cancer models, in particular, cancer invasion models, by proposing a new description of proliferation. More precisely, we extend the classical logistic proliferation function to include a non-local integral term in the growth part. Such modification allows capturing the spatial expansion, i.e. appearance of cells in new locations, without adding artificial terms. To assess the usability of the new model we apply the inverse problem methodology, i.e. we provide a Bayesian inference for unknown parameters and we demonstrate the accuracy of estimators on experimental data on multicellular spheroids growths for three different cells lines [19]. We determine the reliable range of applicability of the newly proposed model. Since the Bayesian approach requires a large number of simulations, we use the fact that, under some assumptions, the proposed model is radially symmetrical, and we transform it into polar coordinates. Finally, we prove the stability of the numerical scheme used to solve the model and we give conclusions.
The structure of the paper is as follows. In Section 2 we introduce the new mathematical model, a non-local proliferation function of the logistic type, and compare it to a previously proposed one.
Then, we present the experimental data that we found best suitable for parameter estimation of simple cellular colony growth. We conclude Section 2 with theoretical consideration on the range of the model applicability. In Section 3, we briefly describe the adopted Bayesian inference methodology and we present the observation model. In Section 4 we prove stability of posterior distributions. Then in Section 5, we show the results of parameter estimation, and finally, in Section 6 we discuss the obtained results and direction of future research. For the interested reader, we include Appendix A and Appendix B, one containing an auxiliary lemma on the Lipschitz continuity of the inverse function, and the second one containing results on the convergence of measure solutions of the considered problem (we also refer to and our companion paper for more analytical details [29]).
Appendix C contains the random walk Metropolis-Hastings algorithm.

Non-local Proliferation Model and Experimental Data
A proper description of cell colony development is a very difficult task. This is mainly because many factors influence those dynamics. Even when considering the increase in cell number alone, in addition to cell division an important role plays access to nutrients and population density. The latter particularly affects the dynamics of larger cell colonies, being significantly less important in the development of colonies at the initial stage. Moreover, the overlap of so many phenomena causes the setting up of an experimental regime suitable for obtaining data for parameter estimation of proliferation function very challenging.
An approach that researchers use frequently to picture the increase in cell number is the logistic function, which is sometimes modified for example by adding volume filling term. However, as already mentioned, the simple local logistic function applied to describe the proliferation of cells within a colony is unable to capture its spatial expansion. This problem is often bypassed by the inclusion of artificial terms, primarily diffusion, and different types of taxis, making the description more phenomenological and more distant from the biological process.
Therefore, we propose a new non-local logistic function to describe the proliferation of cells living within a colony where the integral term is introduced in the growth part to capture the phenomena of the emergence of daughter cells adjacent to proliferating cells, i.e.: where α stands for proliferation rate, and k = k(x) is a kernel function with compact support such that, We assume that k(x) is fixed in time radially symmetric kernel with profile K that is k(x) = K(|x|) for x ∈ R 3 . An interesting issue is the choice of a particular shape of kernel k. In our approach, we choose a kernel corresponding to a normalised characteristic function of a ball, i.e.: (2) where σ k stands for kernel size. Note, that the inhibitory term (1 − n(x, t)) remains local according to the interpretation that the emergence of a cell in a given location is limited by the density of cells in that very place. To the best of our knowledge, a non-local logistic proliferation function given by (1) was not proposed before. However, another type of non-local logistic proliferation function was published previously by Maruvka & Shnerb in 2006, who suggested including an integral term in the inhibitory part [35]. Importantly, in their approach, the colony progression in space is again induced solely by the diffusion term, that they keep in their model.
It is intuitively clear that the proposed model (1) describes the dynamics of the cellular colony whose local maximal density is limited by the carrying capacity of the environment, and spatial progression is driven by the emergence of new cells in the neighbourhood of dividing mother cells.
However, to reliably assess the model and possibly determine its range of applicability we need to refer to real data. Undoubtedly, the best data are those obtained within an experimental regime limiting the influence of phenomena other than proliferation itself to the biggest possible extent.
Despite being quite old, from that perspective probably the best data are the classical data published by Folkman & Hochberg in 1973 showing the evolution of multicellular spheroids in laboratory conditions with the medium being replenished and open space made available [19]. The experiments were carried out for three different cell lines, i.e L-5178Y murine leukaemia cells, V-79 Chinese hamster lung, and B-16 mouse melanoma. Even though the experiments were carried out for three different cell lines the cultivated spheroids experienced the same general growth pattern. Figure 1 shows the dynamics of mean diameter and standard deviation of measured for cultivated spheroids of L-5178Y, V-79, and B-16 cell lines, redrawn from the original paper. At first, the spheroids enlarged exponentially for a few days before the onset of central necrosis and then, for several weeks, continued on linear growth beginning with the appearance of necrosis [19]. After reaching a critical diameter the spheroids experience no further expansion [19]. Although the general growth pattern was the same, the precise age of switching between exponential and linear growths and stabilisation differed for cell lines under investigation. Finally, we would like to mention that in previous literature other interesting approaches are describing the dynamics of multicellular spheroids, for instance, proposed by Byrne and Chaplain [6][7][8] and later analysed in many analytical papers [12,13,20,21].
We now turn to a short theoretical discussion on the possible range of applicability of the model (1). According to Folkman & Hochberg, the exponential growth of spheroids lasts for several initial days only [19]. Then the increase in colony volume is (approximately) proportional to its volume and this is because initially, all cells divide regardless of where they are located in the spheroid.
However, that period in colony development continues as long as the cell number is relatively small and therefore is beyond the scope of applicability of any density model -for that initial period of colony growth, the discrete description is more adequate. Since the appearance of the central necrosis, the proliferation within the spheroid is limited to the outer layer of several rows of cells.
That means that the colony growth becomes (again approximately) proportional to its surface area.
Assuming the range of the kernel the k(x) to be significantly smaller than the radius of the spheroid, then such a scenario is suitable to describe with (1). Finally, to explain the observed cease of growth of spheroids we theorise that, over time, that their dynamics become influenced by processes that are not present at the begging of cultivation, which go beyond the phenomenon of proliferation. This may be, for example, the lysis of the necrotic part of the spheroid or an inhibitory effect resulting from the appearance of catabolites or necrotic debris.
Taking the above into account, we formulate a hypothesis that the model, given by the (1) is suitable to describe the evolution of cell colonies in which the growth is limited to the outer layer of cells. To bolster the hypothesis we confront the model with experimental data presented on Figure 1 and, for the period of quasi-linear growth of the tumour diameter, we perform parameter estimation using the Bayesian inference approach. Although it seems a natural choice, the kernel given by 2 causes some analytical and practical problems. To solve (1) equipped with compactly supported and Lipschitz continuous kernel one could apply an approach based on the numerical scheme called Escalator Boxcar Train (ETB) developed by de Roos [15]. However, the low regularity of the kernel which is not Lipschitz continuous prevents using standard arguments to prove the convergence of the algorithm. Therefore, we use the fact that both initial data and kernel k(x) are spherically symmetrical, and we rewrite (1) using spherical coordinates that gives us (3). We refer the reader interested in the detailed arithmetic of the change of variables to the appropriate section in our companion paper [29].
Let n(x, t) be the solution to (1) with radially symmetric initial condition n 0 (x). Then the radial density p(R, t) defined as p(R, t) = 4πR 2 n((0, 0, R), t) with p 0 (R) = 4πR 2 n 0 ((0, 0, R)) satisfies where the interaction kernel L(R, r) is given by The equation (3) has Lipschitz kernel with only one singularity at zero. Using appropriate weighted norms we have proven that the numerical algorithm based on ETB approach converges in the setting of Radon measures. Since the full proof of this fact goes beyond the scope of this paper again we refer to Appendix B and our analytical paper for more details [29]. Here, let us only remark that the EBT, or more general a particle methods, for (3) boils down to the following. As the solutions to (3) are supported for all R ≥ 0 even if one starts with compactly supported initial conditions, we introduce R 0 > 0 such that p(R, t) is negligibly small for R ≥ R 0 , see Theorems B.1 and B.5 for the precise statement, and we approximate the distribution where x i = i N R 0 and i = 1, ..., N . With these assumptions in place, it is now sufficient to solve the system of ODEs for masses m i (t): where m i (0) are chosen so that An additional benefit of the change of variables is the improvement in computational accuracy as well as the speed up of the numerical simulations, which is particularly important considering that Bayesian inference usually requires thousands of iterations of solutions to the estimated model. We conclude the section with the remark that techniques based on the flat norms on spaces of measures became recently a promising tool for optimal control problems [1,24,39] which may result in the future in combining Bayesian techniques with optimal control.

Observation Model
To estimate parameters of the model (1) we propose to use Bayesian inference. Within this approach, unknown parameters are treated as random variables that can be described with probability distributions. Bayesian inference is based on posterior distribution and the conditional distribution of parameters given the observed data. The posterior distribution, by Bayes theorem, is given by where D denotes collected data and θ denotes a given vector of parameters, whereas Θ is the space of all parameters. To be precise, π(θ|D) is the posterior probability density that is the probability density of θ given data D, π(θ) is the prior probability density, that is the probability density of θ without any knowledge on data, finally, (D|θ) is the likelihood function that quantifies the probability of observing data D, given the parameter θ.
Usually, it is not possible to obtain an analytical formula of the joint posterior distribution π(θ|D) given by (8). A possible way to overcome this difficulty is to use numerical methods, of which the very popular are the Markov chain Monte Carlo (MCMC) methods [3,36]. MCMC methods comprise a whole class of algorithms including one of the most widely used -the Metropolis-Hastings algorithm.
The Metropolis-Hastings algorithm can be used to generate a sample from the posterior distribution π(θ|D), which in turn can be used to estimate the mean of the posterior distribution. The main idea behind that algorithm is to simulate a Markov chain whose stationary distribution is π(θ|D).
This means that for a sufficiently large number of steps, samples from the Markov chain look like the samples form π(θ|D).
The first state θ 0 of the Markov chain is selected according to some chosen "a priori" distribution π(θ). "A priori" distributions reflect our belief about the nature of the estimated parameters. Such a belief may be based on intuition, experience, assumptions, or even a simple guess. Then, the next step in the Metropolis-Hastings algorithm is selecting a candidate for the next state of the Markov chain taking into account the current state θ j . This requires defining a method of sampling the parameter space Θ, i.e. requires defining a probability density g(θ|θ j ), sometimes referred to as the proposal density or jumping distribution, that suggests a candidateθ for the next sample value θ j+1 , given the previous sample value θ j . In the case of an unknown parameter being a number, the probability density g(θ|θ j ) is often chosen to be a normal distribution. Whereas, when unknown parameters are a vector, the probability density g(θ|θ j ) is usually a multivariate normal distribution, which makes the proposing of a candidate for a new state from a current one very simple. Having  Except for the gray rectangles, data are rewritten from Folkman & Hochberg [19].
the candidate for the next state of the Markov chain we have to decide whether to accept it or not. There is no single criterion for doing that, however widely used is the function proposed by Metropolis, i.e.
In the case when g(θ j |θ) = g(θ|θ j ), i.e. the proposal density is symmetrical, (9) simplifies to Using (8) we obtain: Function A, often called the acceptance probability, gives the probability of the candidateθ being accepted as the next state in Markov chain. The Metropolis-Hastings algorithm generates a uniform random number u ∈ [0, 1] and if u ≤ A then sets θ j+1 :=θ, otherwise sets θ j+1 := θ j .
To calibrate the proposed proliferation model we need to set up the link between the data and the theoretical framework namely, we need to define the so-called observation model. To estimate the parameters we use three series of measurements of diameters of multicellular spheroids provided by Folkman & Hochberg [19]. Assuming that all spheroids have almost homogeneous mass, we assume that colony radius at time t determine the sphere containing 95% of the current mass of the spheroid, i.e.
where r(t) denotes the colony radius at time t, whereas ∞ 0 p(r, t) dr stands for colony mass at time t. Continuity of such a defined quantile function strongly depends on measure p. On the subset of measures having density with connected support, it can be shown that the quantile is Lipschitz continuous with respect to the underlying measure. This property is crucial for showing the stability of posterior distribution. In our approach, the solution is given by a discrete measure µ N t , so to make the solution continuous, one may convolve it for instance with Laplace distribution The test comparison computations we conducted indicate that such a regularisation is not needed in practice. The differences between the simulation results are imperceptible and concern distant decimal places. On the other hand, the regularisation requires much more computing power, in particular, it extends the calculation many times, therefore we omit it while performing the parameter estimation.
Despite any efforts, the measurements of the spheroids diameters are of course burdened with an error, therefore we assume that or alternatively where r i o stands for colony radius at measurement performed at time t i , r(t i ) is the actual colony radius at time t i and σ o stands for homogeneous over time measurement error. In conclusion, we perform parameter estimation for three data series where l stands for number of measurements taken into consideration.
The choice of the function describing the initial distribution is to some extent arbitrary. It seems reasonable to assume that the initial function is close to the characteristic function of the ball with the phenomenological modification consisting of mollifying the edge to capture the fact that the cell density on the colony surface is smaller than inside. We assume that it is given by (16) p(r, 0) = 4πr 2 whereσ i and q are chosen so that the radii of the initial colony calculated according to the formula In result we have a vector θ = [α, σ 2 k , σ 2 o , σ 2 i ] of unknown parameters whose coordinates correspond to proliferation rate, kernel size, measurement error, and initial colony radius, respectively. For convenience, and to avoid unnecessary constraints, that are α, σ k , σ o , σ i > 0, we use logarithms of parameters instead of parameters itself. Finally, we need to define the likelihood function (D|θ) whose form follows directly from the assumption (14) and is given by: where r(t i ) is the colony radius at time t i for vector of parameters θ.
To sample the parameter space Θ we choose a random walk Metropolis algorithm, i.e. the proposal distribution is a multivariate normal distribution where s is a step-size, tuned in a way that the acceptance probability of the candidateθ is close to the optimal one [22].
Finally, to complete the description of the model, we need to provide the specific "a priori" distribution π(θ). We assume the time scale of simulations corresponding to the time scale of the experiments conducted by Folkman et al. [19], therefore, we set the time unit to 1 day. For similar reasons, we take 1 mm as the length unit. We assume that "a priori" distributions for all unknown parameters are independent log-normal distributions. When determining the "a priori" distributions of proliferation parameters for particular cell lines, we use the mean doubling time for each cell type, see Table 1 and we assume that cells in a colony proliferate 2.4 times slower. Within this framework, we get the "a priori" proliferation rates equal to 1.4, 1.04, 0.9, for L-5178Y, V-79, and B-16 cell lines, respectively. According to Folkman & Hochberg [19], the proliferating ring is restricted to the outer layer of several cells. Therefore we assume that kernel size is equal to 6 times the cell diameter, which gives us mean values of kernel size equal to 0.06 for L-5178Y and V-79 and 0.09 for B-16 cell lines. We assume that the mean value of "a priori" distributions of initial colony radius is equal to the mean radius of the first measurements, moreover, the false and positive errors are equally probable therefore we set the mean values of measurement errors equal to zero. For α, σ k and σ i we set the standard deviation to 1, whereas for σ o we set it equal to 5 to cover the rather higher uncertainty of observation error over other parameters.

Stability of posterior distribution
Bearing in mind that we approximate the posterior distribution (8) using numerical solutions rather than exact ones, we need to prove that our approximation is indeed close to the actual one. At first, in Subsection 4.1 we propose a regularisation of the quantile function used for colony radius approximation (12) and prove its Lipschitz continuity. We refer the interested reader to the Appendix for the detailed formulation of an auxiliary lemma, which is needed to show that the inverse of the cumulative distribution function satisfies the continuity estimate. Within the Appendix, we also recall the necessary notions from measure theory and we formulate theorems about the existence and uniqueness of measure solutions. In Subsection 4.2 we prove the stability of posterior probability distribution of θ with respect to the EBT approximation of (3).

Regularisation of percentile function and its Lipschitz continuity.
The quantile function (12) used to determine the colony radius at time t is not invertible, which is crucial for further analysis as we retrieve radius from the measure solution to obtain the likelihood function (17).
Hence, we propose the following regularisation. Consider µ ∈ M + (R + ) and continuous function where µ * η is a convolution of the measure µ with the function η, which we call a regularising kernel, and µ T V is the total variation norm defined by (28). We note that such a convolution is a function as well. We impose the additional assumptions on η and initial condition µ 0 that guarantee stability properties of (19).

Lemma 4.1.
Let µ t ∈ M + (R + ) be a measure solution to (3) with initial condition µ 0 and let µ N t be the particle approximation defined in (6)- (7). Then, there is a constant C depending continuously on parameters and the size of initial conditions µ 0 T V such that Proof. Measure solutions to (3) are non-negative and uniformly bounded with respect to · T V , with a constant depending only on the initial condition, time, and parameters, see Theorem B.1 in Appendix. The proof of that theorem can be found in our companion paper [29]. Using (3) we understood in the sense of distributions. Taking convolution with η we deduce Integrating in time we conclude estimates for µ t . To establish estimates for µ N t we observe that (6) implies distributional inequality so that the proof above applies also to µ N t . As µ 0 T V = µ N 0 T V , the proof is concluded. Now, we are in position to prove that on the appropriate set the function F η satisfies Lemma A.1.
For simplicity, we denote byθ = [α, σ k , σ i ]. Note that the constants in the estimates (22) are independent ofθ assuming that α, σ k , σ i are in the certain range of values, that is usually bounded and separated from zero. For the forthcoming consideration, it is convenient to define two sets R and S consisting of solutions to equation (3) and the numerical scheme (6)- (7), respectively. Moreover, to investigate stability properties of radial solutions we use weighted flat norm defined by (30).

Remark 4.4.
The existence of an appropriate R 0 follows from Remark B.3.
Proof. First, we note that ∂ x F η (x, µ t ) = µ t * η ≥ e −CT µ 0 * η ≥ ε κ so that we can always find such a κ uniformly for all elements of R. Existence of such b κ follows from uniform tail estimate (32) in Theorems B.1 and B.5.
Concerning Lemma A.1, the first condition is satisfied. For the second, we write Note that The function z → x −∞ η(y − z) dy is bounded by an L 1 norm of η and Lipschitz continuous as for all z 1 , z 2 we have as η was assumed to be Lipschitz continuous with constant C Lip (η). It follows that and consequently, using Lemma 4.1 we can estimate term A with For term B we observe µ t T V = R + dµ t and ν t T V = R + dν t so that where C comes from Lemma 4.1. By virtue of Young's convolutional inequality we observe that Finally, we note that for all R 0 > 0 we have interpolation inequality For the first term we note that the map [0, R 0 ] r → rψ(r) is bounded with R 0 and Lipschitz continuous with constant (1 + R 0 ). Hence, For the second and third term we use decay estimate (32). Indeed, |r ψ(r)| ≤ 2 and r/2 ≤ e r/2 so Taking supremum over all ψ ∈ BL(R + ) with ψ BL ≤ 1 we conclude the proof of (4.1) which proves where C may depend on ε(µ 0 ), ε(θ), ε(κ), κ, a κ , b κ . Now, as κ < 0.05 we obtain (23) directly from Lemma A.1.

Proof of stability of posterior distribution. Let us remind, that posterior distribution is
given by (8) with likelihood function defined by (17). The actual colony radii r(t i ) are the function of the solution, i.e. r(t i ) = G µt i . Hence, we may write where we added a superscript µ to denote dependence on the measure solution µ. Recall that as in Theorem 4.2, we work in the set R ∪ S of measure solutions obtained with appropriate initial conditions and values of parameters as well as solutions to the numerical scheme.  Proof. Note that log a κ ≤ log(G µt i (θ) ) ≤ log b κ , so the first inequality follows directly from assumptions and formula (24), while the second one follows from the first after noting that Θ π(θ) = 1. Proof. Clearly, the function x → e x for x ≤ 0 is 1-Lipschitz. Moreover The conclusion follows.
Assume additionally that σ o ≥ ε(θ) > 0. Then, there is a constant C such that for all R 0 > 0, Proof. First, we observe that We note that triangle inequality and inequality 0 ≤ e −x ≤ 1 for x ≤ 0 implies Then from Lemma 4.7 we obtain that where a κ and b κ are such that 0 < a κ ≤ G µt i (θ) , G νt i (θ) ≤ b κ . Letting Then, equation (23) implies for a possibly larger constant C. Hence, from (26) we deduce where we applied Lemma 4.5 and Θ π(θ) = 1.
Theorem 4.9 (Stability of posterior distribution with respect to particle approximation). Let µ t ∈ R be a measure solution to (3) with initial condition. Let µ N t = m i (t) δ xi where m i (t) solve system of ODEs (6) and m i (0) are chosen as in (7). Then, if we let .
we have In particular, Proof. We let ν(t) = µ N (t) in (25) and use (33) to conclude the proof.

Simulations Results
To perform Bayesian inference, i.e to predict the growth curve of diameters of multicellular spheroids and to identify parameters related to proliferation rate, kernel size, measurement error, and ini-        growth the MAP estimator seems to be more accurate. Using the MAP estimator we obtain the prediction of diameters dynamics very close to the Bayesian one however, this approach allows us to obtain more accurate parameters for the proposed proliferation function (1). Figure 5 presents the predicted dynamics of diameters for all considered data sets, whose prediction were obtained using the MAP estimator (red curve) and the Bayesian estimator (blue line). Using the MAP estimator we due to the correlation between parameters α and σ k that for the cell line whose linear growth is the slowest becomes more apparent. We speculate, that for such challenging cases might be worth trying more sophisticated algorithms Metropolis-Hastings MCMC, however, the issue goes beyond the scope of the current paper, whose main aim was to propose a new proliferation function suitable to incorporate into models describing solid tumour dynamics. Finally, we observe that estimated parameters in all examples are similar, which bolsters the surmise that the model does not overfit the data. Therefore, our model provides a quite good approximation of reality in the considered time window.

Conclusions
In this paper, we propose a non-local function (1) to describe the proliferation dynamics of cells living within a colony whose growth is restricted to the outer layer of several viable individuals. To estimate the range of applicability of the model, we refer to the experimental data on cancer multicellular spheroids growth provided by Folkman & Hochberg [19]. To deal with the low regularity of the kernel given by (2) as well as to improve solution accuracy and algorithm performance we reformulate the initial model using radial coordinates (3). Then, we performed parameter estimation of the model given by (3)  An interesting question arises, to what extent the proposed description of cell proliferation is suitable to incorporate into more complex cancer models. Whether the introduction of the proposed non-local proliferation function will bring more accurate quantitative predictions of solid tumour growth or not? To answer that question it seems interesting to relate the estimated kernel radii to the distance that oxygen and nutrients can effectively diffuse into living tissue. It is known that the threshold that oxygen can effectively diffuse through tissue is about 0.2 mm [43]. Considering that not only oxygen is needed to keep cells alive, but also nutrients, whose molecules are larger, the distance between capillaries and the necrotic core will be smaller than mentioned 0. we postulate its suitability for describing proliferation in cell colonies whose growth is restricted to outer layers of viable cells. Let us mention that such a scenario is typical for most solid tumours, which, due to the lack of a regular blood vessels network, typical for healthy tissues, develop its blood supply through the process of angiogenesis that results in a pathological capillary network producing numerous necrotic regions.
While analysing proliferation parameters for considered cell lines obtained with the MAP estimator becomes conspicuous that the value 1.7264 obtained for the L-5178Y cell line is noticeably larger than the values obtained for V-79 and B16 that are 0.3603, and 0.3616, respectively. The fact becomes more comprehensible if one considers also the dynamics of the entire colonies. Spheroids composed of L-5178Y cells grow much faster and reach a diameter of about 4 mm after only 30 days.
The growth is very fast but lasts shortly. What's more, the estimated value corresponds to the cell doubling time of about 10 hours, which is exactly the value reported in databases [32]. In summary, it seems that in the initial growth of the L-5178Y spheroid, the presence of neighbourhood cells Summing up our work we state that comprehensive calibration of complex models describing the dynamics of the cancer disease in vivo seems for the moment to be out of reach, first of all, due to the lack of relevant data but also due to computational complexity of such tasks. Therefore, it seems appropriate to create at first partial, properly calibrated models, that describe phenomena contributing to cancerogenesis and then combine them into more complex models to get more quantitative insight into the pathology of cancer development. Appendix A. Lipschitz continuity of the inverse function A simple lemma, which is in fact an inverse function theorem with a parameter. We use this lemma while proving that the inverse of the cumulative distribution function satisfies the continuity on (a, b). Note that we are taking advantage of the fact that the problem under consideration is now one-dimensional since for multidimensional cases only local invertibility is true.
Then, for all y ∈ (a, b), Proof. Note that standard inverse function theorem implies that for all s ∈ S we have (f −1 s ) ∞ ≤ 1 δ . Hence, we can estimate: Appendix B. Theory of measure solutions and convergence of particle method Using the theory of measure approach for biological problems is becoming increasingly popular as it is very intuitive and convenient for both analytical and practical numerical simulation viewpoints.
In general, measures assign real values to all measurable subsets of the considered space -we say that those real values are measures of sets. Naturally, non-negative measures are useful for describing various observable quantities, such as age, distribution, or density. Importantly, the spaces of measures are vector spaces, which means in particular that the difference of two measures is again a measure from the same space. The numerical approximation we use to simulate the solutions to the model (3) is based on distributions not necessarily having densities with respect to the Lebesgue measure. Therefore, we reformulate our model for a generalised class of solutions in the space of non-negative Radon measures, cf. [9,17,18,[26][27][28]42]. Within the new approach µ t (A) is a measure such that, for every measurable set A ∈ R + , µ t (A) is the mass of cells being at time t at a distance from the centre of the coordinate system belonging to the set A, i.e. (27) µ t (A) = A p (r, t) dr.
Using the ETB method we approximate the solution to (3) assuming that each cohort represented by a Dirac mass corresponds to a mass of cells belonging to the set A, see (5).
We denote by M(R + ) the space of all bounded and signed Radon measures on R + whereas M + (R + ) stands for its subset, consisting of non-negative measures. Let's note that for µ ∈ M(R + ) we have unique Hahn-Jordan decomposition µ = µ + − µ − where µ + , µ − ∈ M(R + ). To work in the spaces of measures one needs the notion of a norm. The total variation norm is defined by which can be thought of as the total mass of µ. Moreover, we denote with · BL * the flat norm where space of bounded Lipschitz functions BL(R + ) is given by where C is a constant depending continuously on parameters 0 < α, σ k , σ i . Moreover, we have the following decay estimate: if µ 0 satisfy R + e x dµ 0 (x) < ∞ then In particular,  (32) allows assuming that the support of µ t is finite.

Remark B.4.
It is classical to apply flat norm in problems related to the convergence of particle method cf. [10,11,[25][26][27]. However, our kernel L(R, r) defined by (4) is singular at R = 0 or r = 0 and therefore it does not satisfy the assumptions of the previous works. Despite this difficulty in our theoretical paper, we were able to prove convergence in the weighted flat norm [29].
We also recall the result on the convergence of the particle method and certain properties of particle approximation, see [ where C is a constant depending continuously on parameters 0 < α, σ k , σ i .