Bayesian analysis of the effect of exosomes in a mouse xenograft model of chronic myeloid leukemia

1.   Introduction

The National Cancer Institute (NCI) defines cancer as a disease in which the cells of an organism grow uncontrollably and at an accelerated rate, altering cellular homeostasis ^[1]. The study of cancer involves monitoring the growth of tumors, usually using mathematical and/or statistical tools. In the first decade of the 20th century the use of systems theory to study biological processes became popular; von Bertalanffy was a pioneer in this area, he was the author of the book "Systems Theory." Published in 1968, he proposed a growth equation that bears his name ^[2]. Mathematical modeling has gained popularity in biological processes and in particular the growth of tumors.

Some of the classic tools for modeling tumors are the growth equations, such as exponential growth proposed by Thomas Malthus in 1798, Gompertzian growth proposed by Benjamin Gompertz in 1825, logistic growth proposed by Verhulst in 1838, among others ^[3]. More sophisticated models have been proposed throughout history. In the 1990s, James Murray studied the effect of chemotherapies on the spread of tumors by means of computed tomography ^[4]; at the end of the 20th century, Chaplain proposed models for tumor vascularization ^[5]. Currently, work is being developed by using artificial intelligence for the simulation of tumor microenvironments ^[6], among other works that add to the battle against cancer, which seems to be eternal.

Chronic myeloid leukemia (CML) is a hematological neoplasm of slow progression that originates in the bone marrow, caused by the Philadelphia (Ph) chromosome ^[7]. The Ph chromosome is an abnormally short chromosome 22 that is one of the two chromosomes involved in a translocation with chromosome 9, producing an overproliferation in the myeloid cell line, and generating imbalance in the cell populations of the blood. The symptoms of the latter are the following: tiredness without reason, low body weight, pallor, and skin spots. Most patients with CML do not present symptoms until the most advanced stages of this cancer. Lin et al. reported an annual incidence of 34, 179 cases and 24, 054 deaths related to CML worldwide in 2017 ^[8]. By 2023, the American Cancer Society estimates that one in 526 people will suffer in their lifetime from CML in the United States ^[9]. This disease is most often diagnosed in people over 60 years of age with higher percentage in men.

One treatment available for CML, that decreases the rate of cancer cell proliferation, is Imatinib (IM), which is an orally administered drug used to block the enzyme tyrosine kinase involved in cancer cell proliferation, thus slowing cancer progression. Retrospective studies by Adattini et al. ^[10] on 86 patients with CML medicated with IM in the period 2001–2018 showed that treatment with IM induced a positive effect on patient survival. Brümmendorf et al. ^[11], reported an estimated 5-year overall survival rate greater than 90%. IM treatments have been extensively investigated and are considered safe. An in vivo study in a mouse xenograft model showed that exosomes promoted the proliferation and decreased the sensitivity of CML cells to IM, resulting in treatment resistance ^[12].

The CML originates in the bone marrow, in the presence of human bone marrow mesenchymal stem cells (hBMMSC), these cells release nanovesicles of 30–150 nm in diameter called exosomes (hBMMSC-Exo), endowed with a diverse load of proteins, lipids and nucleic acids ^[13]. hBMMSC-Exo can stimulate tumor initiation and therapy resistance in cancer cells ^[14], due to their cell-cell communicative function, i.e., the rate of cancer cell proliferation increases in the presence of exosomes from hBMMSC. For this reason, it is necessary to explore how the proliferation dynamics develop in the presence of hBMMSC-Exo and in IM-mediated therapy.

A classic way of quantifying cancer cell proliferation is through the volumetric tumor growth, for which mathematical models based on already established differential equations can be used. These models describe the number of cells or tissue volume as a function of time, and they have the following form:

where $f$ is a function that characterizes the growth dynamics, $V(t)$ is the tumour volume at time $t$, and $\boldsymbol{\theta}$ is a vector of parameters. Logistic and Gompertz models are widely used. These are expressed mathematically in Table 1.

In these cases $\boldsymbol{\theta} = (r, K)$, the parameters $r$ and $K$ represent the growth rate in cubic millimeters per unit of time ($\text{mm}^3 t^{-1}$) and the carrying capacity in $\text{mm}^3$, respectively.

The parameter vector $\boldsymbol{\theta}$ can be modified with the factors $(1+e_j)$ or $(1-i_j)$ with $j = 1, 2$, when estimating the effect of a growth stimulant or inhibitor, such as Ex and IM, respectively. $e, i \in [0, 1)$ are interpreted as percentages of the stimulant or inhibitory effect on the parameters $r$ and $K$. Therefore, the growth models in Table 1 can be rewritten with effects, as shown in Table 2.

The aim of this work is to quantitatively estimate the effects of exosomes (hBMMSC-Exo) and Imatinib (IM) on the growth rate and carrying capacity of CML. Five specific objectives are proposed: First, to estimate the parameters of the CML tumor growth; second, to estimate the parameters of the hBMMSC-Exo effect in CML; third, to estimate the parameters of the IM effect in CML; fourth, to estimate the parameters of the combined Ex-IM effect in CML; and fifth, to determine the model with the best predictive performance in the CML observations.

2.   Materials and methods

2.1.   Experiment design: Tumor growth with treatments

To perform parameter estimation of the effect of hBMMSC-Exo and Imatinib on CML, we used the observations made by Zhang et al. ^[12] with the human CML cell line K562. Four tumor growth assays were performed by using twenty female BALB/c-un mice that were divided into four groups. Each group of five mice received a subcutaneous injection of $200$ $\mathrm{ml}$ into the front of the right rear of the mice containing K562 cells, the following four treatments were performed for 40 days. Tumor volume was measured every four days for a total of ten measurements per treatment. Tumor growth was evaluated by tumor volume, which was calculated by using the modified ellipsoidal formula: $V = 0.5 \times length\times width^2$. The detailed procedure of the xenograft protocol can be found in Zhang et al. ^[12]. The data are recorded in Figure 1.

2.2.   Bayesian estimation

2.2.1.   Statistical model

Let $\mathbf{y} = (y_1(t_1)$, $y_2(t_2), \dots, y_n (t_n))$ be the vector of $n$ observations of tumor volume over time. The following statistical model was considered:

where:

● $y_{i}(t_i)$ is the dependent variable of the model and represents the tumor volume at time $t_i$. for the four treatments of the study.

● $\boldsymbol{\theta}$ is the vector of parameters corresponding to each stage of the Bayesian estimation.

● $G(\mathcal{M}(t_i, \boldsymbol{\theta}))$ is the solution of the model $\mathcal{M}$. For practical purposes, we obtained a numerical solution using the Runge-Kutta fourth order method.

● $\varepsilon(t_i)$ is the random error at time $t_i$, the errors are independent for each time, normally distributed with zero mean and variance $\sigma^2$.

2.2.2.   Likelihood function

Considering the Gaussian assumption $\varepsilon(t_i) \sim N(0, \sigma^2)$, then $y_{i}(t_i) \sim N\left(G(\mathcal{M}(t_i, \boldsymbol{\theta})), \sigma^2\right)$. For such a reason, the likelihood function of $\boldsymbol{\theta}$ based on $\mathbf{y}$ is given by the following equation:

2.2.3.   Posterior distribution

Bayesian statistics allows incorporating information or knowledge of the parameters to the inference process. This information is specified by the researcher by means of a prior distribution denoted by $P(\boldsymbol{\theta})$. The prior distribution takes into account the range of values that the parameter takes and can assign a higher probability to a subset of those values or, in the absence of information, assign the same probability to all values. For our case, the tumor growth rates $r$ should be in the interval of $[0, 1]$ and $K \geq 0$, and do not exceed $1500 \text{ mm}^3$ over the 40 days of the study.

The Bayesian inference is based on the so-called posterior probability distribution, which is defined by Bayes' theorem as follows:

where

● $P(\boldsymbol{\theta} \mid \mathbf{y})$ is the posterior probability of $\boldsymbol{\theta}$ given a set of observations $\mathbf{y}$.

● $L(\boldsymbol{\theta} \mid \mathbf{y})$ is the probability of the observations $\mathbf{y}$ with specific value of the vector $\boldsymbol{\theta}$, that is, it is the likelihood.

● $P(\boldsymbol{\theta})$ is the prior distribution.

● $P(\mathbf{y})$ is a normalizing constant defined as $\int_p L(\boldsymbol{\theta} \mid \mathbf{y}) P(\boldsymbol{\theta}) d \boldsymbol{\theta}$.

Note that the expression (2.3) is well defined if $P(\mathbf{y})\neq0$, $P(\mathbf{y})$ is constant, therefore, we can rewrite it as

The Bayesian inference is based on the posterior distribution that summarizes the information in the data and the prior information about the parameters. The mean or median of the posterior distribution is very important to Bayesian inference, and by definition, they are obtained as a multidimensional integral of the posterior distribution

however, most of the time this integral is not possible to obtain analytically. This problem has been solved using numerical methods such as Monte Carlo integration. A better approach is to use Markov chains Monte Carlo (MCMC) that give the name to a group of algorithms to get a sample from a probability distribution using a Markov chain. A Markov chain is a discrete stochastic process in which the probability of an event depends only on the previous event and has the posterior distribution as its equilibrium distribution. Some examples of MCMC are the Metropolis-Hastings algorithm and the Gibbs sampling. A more sophisticated method and with better performance is the Hamiltonian Monte Carlo (HMC).

2.2.4.   Hamiltonian Monte Carlo (HMC)

We use the Hamiltonian Monte Carlo (HMC) sampler. HMC has proven to be a more efficient sampler than the traditional ones. Its acceptance rate is a little less than twice the acceptance rate of the Metropolis-Hastings algorithm ^[15]. This sampling technique is based on Hamiltonian mechanical physics to sample high dimensionality distributions. The reader is referred to Betancourt et al. ^[16] for a detailed description of the HMC technique.

2.2.5.   Chain convergence criteria

Most MCMC algorithms, including HMC, have the desirable property to have the posterior distribution as its equilibrium distribution. The actual problem is to determine how many iterations are needed to converge within an acceptable error. This is the reason why we need convergence diagnostics to know if we have enough iterations to have a good approximation to the posterior distribution that will be used for Bayesian inference.

To favor convergence we use four chains with a random init, and the following convergence diagnostic criteria:

● Traceplot and density compose an empirical approach for convergence control consisting of graphics of the simulated chain output to detect non-stationary behaviors ^[17]. The traceplot is a time series plot of the iteration number vs the realizations of the Markov chain at each iteration. There is convergence when the graph shows good mixing across chains. The density plot is a non-parametric estimate of the density of each chain. There is convergence when the densities are similar.

● Numerical diagnostic Gelman-Rubin's diagnostic compares two estimations of the variance, inter-chain and intra-chain. The ratio of the estimation of the variance gives two statistics: the potential scale reduction factor (PSRF) (also called Rhat, $\widehat{R}$) and its credible interval. There is convergence of the chains when the PSRF and the upper limit of its credible interval are close to one and less than 1.2, respectively ^[18].

2.2.6.   Bayes' estimator

To make inferences, the following Bayes estimator is used

with the quadratic loss function $L(T, \boldsymbol{\theta}) = (T- \boldsymbol{\theta})^2$, we found the minimum at $T = E\left[\boldsymbol{\theta} \right]$, i.e., the Bayesian estimator obtained is the posterior mean $T^B = E(\boldsymbol{\theta} \mid \mathbf{y})$. We used this estimator to make inferences on the posterior densities. The 95% credible interval (CrI) was calculated using the 2.5th and 97.5th percentiles.

2.2.7.   Model comparison

Leave-one-out cross-validation (LOOCV) is a technique used to evaluate the predictive performance of models, and it is mainly used for model comparison. LOOCV is the extreme case of k-fold cross-validation when $k = n$ data partitions are used. The data used only has $n = 9$ observations, that is, it is a small sample. Using $k < n$ partitions implies eliminating at least two observations, which leaves an even smaller sample that can produce unreliable model parameter estimates. This technique, in a Bayesian setting, fits the model with the exclusion of one observation and estimates the expected log pointwise predictive density (ELPD) for the data out of the sample, in practice for a set $y_1, \dots, y_n$ of independent data of $\boldsymbol{\theta}$. The ELPD estimate is defined by the following expression:

The reader is referred to Vehtari et al. ^[19] for a detailed description of the Bayesian LOOCV technique. In practical terms the model with the highest $\text{elpd}_{\text{loo}}$ value is the model with the highest predictive performance.

2.3.   Statistical software

This work was done in the Julia programming language ^[20] using the libraries DifferentialEquations.jl ^[21] and turing.jl ^[22]. Parameter estimation of each model was done using four chains of 3500 iterations of the HMC, and a burning of the first 1000 samples.

The convergence diagnostics, Gelman-Rubin and trace plots, were performed by default in the turing.jl package.

The posterior distributions were fitted by using Julia's Distributions.jl package ^[23].

Once the models were fitted, Bayesian LOOCV was calculated using Julia's Arviz.jl package ^[24] which performs an exploratory analysis of Bayesian models.

2.4.   Bayesian estimation stages and prior distributions

The Bayesian estimation of the parameters of a model for each data set comprises four stages.

1) First stage: The parameter vector $\boldsymbol{\theta} = (r, K)$ for K562 tumor growth is estimated by using the models in Table 1. Non-informative prior distributions are proposed in A1, see Table 3.

2) Second stage: We estimated the parameter vector $\boldsymbol{\theta} = (r, K, i_1, i_2)$ for K562IM tumor growth by using the models with inhibitory effect from Table 2. The prior distributions for $r$ and $K$ in the second stage use the posterior distribution obtained in stage one. The posterior distribution for $r$ is a Beta distribution that is a generalization of the uniform prior from stage one, with hyperparameters obtained by fitting the Beta distribution to the posterior distribution of $r$ in stage one. The posterior distribution for $K$ is a normal distribution because the posterior distribution obtained in stage one has a symmetric shape and the hyperparameters are obtained by fitting a normal distribution to the posterior distribution of $K$ in stage one. Furthermore, since $K$ is the carrying capacity, the possible values of the normal distribution are restricted to the interval $(10, 2500)$ because 10 is a value close to the initial volume, and 2500 is twice the maximum observed volume. For the parameters $i_1$ and $i_2$ we proposed non-informative priors in A2 that correspond to a uniform distribution in the interval $(0, 1)$; see Table 4.

3) Third stage: We estimated the parameter vector $\boldsymbol{\theta} = (r, K, e_1, e_2)$ for K562EX tumor growth by using the models with stimulating effect from Table 2. Prior distributions for $r$ and $K$ were selected in a similar way to the second stage. For the parameters $e_1$ and $e_2$ we proposed non-informative priors in A3 that correspond to a uniform distribution in the interval $(0, 1)$; see Table 4.

4) Fourth stage: Identical to the second stage, but with the K562IMEX observations.

An overview of the methodology is presented in Figure 2.

Remark: Each stage of the estimation requires a prior distribution for $\sigma$. In all cases it will be Inverse-Gamma (0.001, 0.001). The prior distributions of the first stage can be observed in Table 3.

3.   Results

3.1.   First stage

We show the fitting of Gompertz and logistic models to K562 tumorigenic cell growth (see Figure 3). Bayesian estimates of parameters $r$ and $K$ are presented in Table 5. The HMC converged as indicated by traces, densities, and the Gelman-Rubin statistic, which was less than 1.2 ^[18] for all the parameters in the models (Figures A2 and A2).

Positive $r$ values for the Gompertz model and logistic model indicate a tumor growth of $0.06$ and $0.19$ volume units per time unit. The $K$ values represent the mean value of the maximum volume limit at the 40th day. Credibility intervals at a level of 95% for $K$ of the models do not overlap, indicating a significant difference in the estimated carrying capacity between the models.

Model selection: The Gompertz model obtained a value of $\text{elpd}_{\text{loo}} = -59.34$ and the logistic model a value of $\text{elpd}_{\text{loo}} = -59.99$. It can be stated that the Gompertz model presents a higher predictive performance for the dynamics of K562 tumor growth, however, both models present a statistically equal predictive performance.

In Figure 3, the solid lines closely follow the observed data points. Although the Gompertz model performs slightly better in terms of predictive accuracy, the difference is small. This suggests that both models reasonably capture the overall trend and underlying processes of tumor growth. The values of $r$ and $K$ obtained indicate that the logistic model suggests faster growth but a lower carrying capacity for the K562 tumor volume compared to the Gompertz model, which predicts slower growth but a higher carrying capacity. The posterior distributions of $r$ and $K$ will be used as a priori distributions in stages 2, 3 and 4 to estimate the effects of the treatments.

3.2.   Second stage

We show the fit of Gompertz and logistic models to K562 tumorigenic cell growth with IM (see Figure 4). Bayesian estimates of parameters $r$, $K$, $i_1$ and $i_2$ are presented in Table 6. The HMC converged as indicated by traces, densities and the Gelman-Rubin statistic, which was less than 1.2 ^[18] for all the parameters in the models (Figures A3 and A4).

The estimates of $r$ and $K$ for the models with inhibition remain significantly the same as the estimates obtained in the first stage. The dynamics of tumor growth are affected by drug treatment and the effects are represented in parameters $i_1$ and $i_2$ for each model.

Model selection: The Gompertz model with inhibition obtained a value of $\text{elpd}_{\text{loo}} = -41.37$ and the logistic model with inhibition obtained a value of $\text{elpd}_{\text{loo}} = -34.62$. It can be stated that the logistic model presents a higher predictive performance for K562IM tumor growth dynamics.

For tumor growth dynamics under the effect of Imatinib drug (K562IM), the logistic model with inhibitor presents higher predictive performance, the drug generates an inhibitory effect of $i_1 = 0.06$ on the growth rate of normal tumor volume $r = 0.19$, which means that in the presence of the drug, the growth rate is 93.5% slower. The carrying capacity $K = 949.22$ for the logistic model is affected by $i_2 = 0.90$, which indicates that the Imatinib treatment reduces the tumor volume to approximately 10% in the first 40 days.

3.3.   Third stage

We show the fit of Gompertz and logistic models to K562 tumorigenic cell growth with exosome (see Figure 5). Bayesian estimates of parameters $r$, $K$, $e_1$ and $e_2$ are presented in Table 7. The HMC converged as indicated by traces, densities, and the Gelman-Rubin statistic which was less than 1.2 ^[18] for all the parameters in the models (Figures A5 and A6).

The estimates of $r$ and $K$ for the models with stimulation remain significantly the same as the estimates obtained in the first stage. The dynamics of tumor growth are affected by exosomes treatment. The effect are represented in the parameters $e_1$ and $e_2$ for each model.

Model selection: The Gompertz model with stimulus obtained a value of $\text{elpd}_{\text{loo}} = -64.82$ and the logistic model with stimulus obtained a value of $\text{elpd}_{\text{loo}} = -58.45$. It can be stated that the logistic model with stimulus presents a higher predictive performance for K562EX tumor growth dynamics.

Tumor growth under the effect of hBMMSC-Exo exosomes (K562EX), the logistic model with stimulation presents higher predictive performance against the Gompertz model with stimulation. The presence of hBMMSC-Exo generates an effect of $e_1 = 0.08$ on the original rate $r = 0.1975$, which means that the growth rate of tumor volume is 0.83 faster. The effect $e_2 = 0.40$ on the original carrying capacity $K = 948.08$ indicated in the presence of hBMMSC-Exo the tumor volume is expanded to 40% in the first 40 days, which allows us to conclude that the growth in tumor volume will be faster and potentially greater than the tumor growth without exosomes.

3.4.   Fourth stage

We show the fit of Gompertz and logistic models to K562 tumorigenic cell growth with IM and exosomes (see Figure 6). Bayesian estimates of parameters $r$, $K$, $i_1$ and $i_2$ are presented in Table 8. The HMC converged as indicated by traces, densities and the Gelman-Rubin statistic, which was less than 1.2 ^[18] for all the parameters in the models (Figures A7 and A8).

The estimates of $r$ and $K$ for the models with stimulus remain significantly the same as the estimates obtained in the first stage. The dynamics of tumor growth are affected by the combined drug treatment and exosomes. The effect is represented in the parameters $i_1$ and $i_2$ for each model.

Model selection: The Gompertz model obtained a value of $\text{elpd}_{\text{loo}} = -63.14$ and the logistic model obtained value of $\text{elpd}_{\text{loo}} = -59.00$. It can be stated that the logistic model presents a higher predictive performance for the dynamics of K562 tumor growth and K562 cell growth with IM and exosomes. Once again, the difference between the values of $\text{elpd}_{\text{loo}}$ is not much and both models are good in their predictive performance.

For Imatinib-exosomes mixture treatment (K562IMEX), the logistic model with inhibitor presents higher predictive performance. The mixture generates an inhibitory effect of $i_1 = 0.17$ on the growth rate of tumor volume $r = 0.19$, which means that in the presence of the mixture treatment, the growth rate is 0.82 slower. The carrying capacity $K = 984.92$ for the logistic model are affected at $i_2 = 0.045$. The Imatinib-exosome mixture treatment reduced tumor volume to 95%.

4.   Conclusions

The Bayesian approach used provides the advantage of working with few observations of tumor growth and the possibility of using prior distributions plays a crucial role in the estimation of the effects generated by the drug and exosomes.

As expected, the Imatinib treatment reduces the tumor volume but the presence of exosomes (hBMMSC-Exo) in the tumor microenvironment leads to a significant increase in the growth of chronic myeloid leukemia xenografts in mice. Furthermore, the exosomes partially reduce the effectiveness of the Imatinib drug in tumor growth reduction. These quantitative results align with the qualitative observations made by Zhang et al. ^[12]. Exosomes found in the tumor microenvironment appear to play a role in facilitating intercellular communication of cancer-related signals.

The conclusions obtained are limited to the xenograft protocol studied. It is worth mentioning that the study of tumors is a complicated phenomenon involving biochemical signals, among other processes that are not addressed in this work.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work has been partially supported by the Consejo Nacional de Humanidades, Ciencias y Tecnologías (CONAHCYT)-Ciencia de Frontera "Ecuaciones de evolución, sus estados estacionarios y su comportamiento asintótico con aplicaciones en Física y Biología" (project N. 684340). Martínez-Fonseca is indebted to CONAHCYT for fellowship that enabled him to pursue graduate studies for the degree of Maestría en Matemáticas Aplicadas.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

In this Appendix, the convergence criteria of the sigmoidal growth models are discussed. The chains of model parameters show good mixing and the $\widehat{R}$ values indicate the convergence of the chains, as shown in Figures A1 to A8.

A1. Markov chains and parameter density for K562 tumorigenic cell growth

A2. Markov chains and parameter density for K562 tumorigenic cell growth with IM

A3. Markov chains and parameter density for K562 tumorigenic cell growth with exosomes

A4. Markov chains and parameter density for K562 tumorigenic cell growth with exosomes