Improving personalized tumor growth predictions using a Bayesian combination of mechanistic modeling and machine learning.

Background: In clinical practice, a plethora of medical examinations are conducted to assess the state of a patient's pathology producing a variety of clinical data. However, investigation of these data faces two major challenges. Firstly, we lack the knowledge of the mechanisms involved in regulating these data variables, and secondly, data collection is sparse in time since it relies on patient's clinical presentation. The former limits the predictive accuracy of clinical outcomes for any mechanistic model. The latter restrains any machine learning algorithm to accurately infer the corresponding disease dynamics.
Methods: Here, we propose a novel method, based on the Bayesian coupling of mathematical modeling and machine learning, aiming at improving individualized predictions by addressing the aforementioned challenges.
Results: We evaluate the proposed method on a synthetic dataset for brain tumor growth and analyze its performance in predicting two relevant clinical outputs. The method results in improved predictions in almost all simulated patients, especially for those with a late clinical presentation (>95% patients show improvements compared to standard mathematical modeling). In addition, we test the methodology in two additional settings dealing with real patient cohorts. In both cases, namely cancer growth in chronic lymphocytic leukemia and ovarian cancer, predictions show excellent agreement with reported clinical outcomes (around 60% reduction of mean squared error). Conclusions: We show that the combination of machine learning and mathematical modeling approaches can lead to accurate predictions of clinical outputs in the context of data sparsity and limited knowledge of disease mechanisms.

Introduction

Advances in patient care have led to the availability of large amounts of data, generated by typical examinations, such as blood sample analysis, clinical imaging (e.g., CT, MRI), and biopsy sampling, as well as by innovative ‘-omics’ sequencing techniques[1,2]. Such clinical data are the cornerstone in the practice of personalized medicine and specifically in the field of oncology[3,4]. However, this abundance of information comes with multiple issues related to data exploitation and synthesis towards the prediction of pathology dynamics. In particular, we identify the following two major challenges: (C1) First, knowledge of the regulatory mechanisms underlying clinical data is largely lacking, and (C2) second, patient data collection is usually sparse in time, since patient clinical visits/examinations are a limiting factor.

Regarding the challenge (C1), scientists have been long supported by the use of mathematical modeling as a tool to identify causal relationships in the experimental and clinical data, particularly in cancer treatment [5-7]. Mathematical models allow to propose and test biological hypotheses, analyze the sensitivity of observables with respect to biological parameters, and provide insights into the mechanistic details governing the phenomenon of interest[8-10]. Although these models can be extremely powerful both in predicting system responses and suggesting new experimental directions, they require adequate knowledge of the underlying biological mechanisms of the analyzed system. Typically, this knowledge is not complete, and only for a limited portion of the involved variables the corresponding mechanistic interactions are sufficiently known. Therefore, even though mathematical models provide a good description of a simplified version of the associated system dynamics, they do not always allow for accurate and quantitative predictions.

On the other hand, machine learning techniques are suitable to deal with the inherent complexity of biomedical problems, but without caring for the knowledge of the underlying interactions[11]. While mathematical models rely on causality, statistical learning methods identify correlations among data[12]. This approach allows to systemically process large amounts of data and infer hidden patterns in biological systems. As a consequence, machine learning-based techniques can provide valuable predictive accuracy upon sufficient training, but do not typically allow for any mechanistic insight into the investigated problem[13]. The overall understanding of the fundamental system dynamics becomes almost impossible, as the chance to generalize the ‘learnt’ system behavior. The latter issue is further exacerbated by the (C2) challenge that has to be faced, related to the sparseness of clinical data. In particular for a single patient, such information is only available at a few time-points, corresponding to clinical presentation. To face the two mentioned challenges with the final aim of improving personalized predictions, we propose a novel—to the best of our knowledge—Bayesian method that combines mathematical modeling and statistical learning (BaM3). As a proof-of-concept, the proposed method is tested on a synthetic dataset of brain tumor growth. We analyze the performance of the new approach in predicting two relevant clinical outcomes, namely tumor burden and infiltration. When comparing predictions from the mechanistic model with those from the BaM3 method, we obtain improved predictions for the vast majority of virtual patients. We also apply the approach to a clinical dataset of patients suffering from chronic lymphocytic leukemia (CLL). The BaM3 method shows excellent agreement between the predicted clinical output and the reported data. Finally, as an additional test case, we show how the proposed methodology can be used to assess the time-to-relapse (TtR) in a dataset of ovarian cancer patients.

Methods

Formal definition of BaM3

We start by assuming a random variable (r.v.) triplet (Y, X, X) that denotes the system’s modelable X, unmodelable X variables/data (e.g., patient’s age or sex, results of different ‘-omics’ techniques, etc.) and the associated observed clinical outputs Y. We then introduce t0 as the clinical presentation time of a patient at which the patient-specific r.v. realizations are obtained. The overall goal of the method is to predict the patient’s clinical outputs by an estimate at a certain prediction time t. The true clinical outputs of the patient will be denoted as y. Moreover, we consider the existence of an N-patient ensemble dataset (y, x, x). In this dataset all the variables (i.e., modelables, unmodelables, and clinical outputs) are recorded at the time of diagnosis t, which might differ from one patient to another. Both t0 and t are calculated from the onset of the disease. We introduce two distinct times to account for the variability of the disease stage among different patients (t) and the time at which a specific patient is presented to the clinic (t0) (see the corresponding Fig. S1).

The core idea of the method is to consider the predictions of the mathematical model as an informative Bayesian prior of the posterior distribution . We can prove that:The implementation of the BaM3 method therefore reduces to the calculation of the aforementioned probability distributions. Although the prediction of the probability distribution function (pdf) of the clinical outputs is rather straightforward for the mathematical model, obtaining the pdf of the patient’s unmodelable data is not trivial. To retrieve the latter, we use a density estimator method upon the patient ensemble dataset to derive p(Y, X), and then consider the patient-specific realization . For further details about method derivation and estimators of performance, see Supplementary Note 1.

Testing the method on synthetic glioma growth

The equations of the selected mathematical model[14,15] (‘full model’) describe the spatio-temporal dynamics of tumor cell density (c), oxygen concentration (n), and vascular density (v) in the context of glioma tumor growth. The full model includes the variation of cell motility and proliferation due to phenotypic plasticity of tumor cells induced by microenvironmental hypoxia [16-20]. It also accounts for oxygen consumption by tumor cells, formation of new vessels due to tumor angiogenesis, and vaso-occlusion by compression from tumor cells[15,21,22]. We generate N = 500 virtual patients by sampling the parameters of the full model from a uniform distribution over the available experimental range. We consider the tumor cell spatial density c to be the modelable variable. Moreover, we treat the integral over the tissue of oxygen concentration and vascular density, denoted as , respectively, as the unmodelable quantities. Starting from the same initial conditions, we simulate the behavior of each virtual patient for 3 years, storing the values of all variables at each month. As sketched in Fig. 1, we use the modelable variable to setup a mathematical model. In particular, we take c(x, t0) at a specific time-point, the clinical presentation time t0, and use it as the initial condition for a Fisher-Kolmogorov equation[23-27] (‘FK model’). We use this model to predict tumor behavior at a specific time in the future, the prediction time t. For each simulated patient we calculate the tumor size (TS) and infiltration width (IW). In parallel, for each patient we evaluate the diagnosis time t as a random number in the interval [t0 − 6, t0 + 6] (in the unit of months), and collect the values of modelables, unmodelables, and clinical outputs at this time to build the patient ensemble. Given the patient-specific modelable and unmodelable variables (c(x, t0) and , respectively) at the clinical presentation time t0, the BaM3 method therefore produces the probability of observing the TS and IW at a specific prediction time t.

Fig. 1

Workflow for generating the synthetic data.

The full model is initialized at t = 0 and used to simulate the spatio-temporal variation of tumor density c, oxygen concentration n, and functional tumor-associated vasculature v. For each virtual patient, two clinical outputs are tracked, i.e., the infiltration width (IW) and tumor size (TS), and two unmodelables recorded, i.e., the oxygen and vasculature integral over the tissue and , respectively. These quantities are generated at the time of clinical presentation t0 every time the method is applied. To generate the patient ensemble, they are also generated at the diagnosis time t. The latter is assumed to be a random time between t0 ± 6 months. Then, the spatial profile of tumor concentration at time t0 is used as the initial condition for the Fisher-Kolmogorov (FK) model, which is in turn used to simulate tumor growth until the prediction time t.

Mathematical models for glioma growth

The system variables are the density of glioma cells c(x, t), the concentration of oxygen n(x, t), and the density of functional vasculature v(x, t)[14,15]. For simplicity we consider a one-dimensional computational domain. We normalize the system variables to their carrying capacity and write the system asHere is a sigmoidal function ((1 + exp(b*(x − x0)), with b > 0 being a constant) allowing for tumor angiogenesis in hypoxic conditions, i.e., for n < n0 where n0 is the hypoxic oxygen threshold. Then, the functions α = α(n) and β = β(n) account for the dependence of cellular motility and proliferation on the oxygen level, respectively[16,17,19]. They are defined as:When the oxygen level is fixed to the maximum level n = 1 in the tissue α = α0 and β = β0, so that the equation for c reduces towhich we denote in the rest of the manuscript as the Fisher-Kolmogorov (FK) model for tumor cell density. We remark that Eq. (7) has been extensively used to predict untreated glioma kinetics based on patient-specific parameters from standard medical imaging procedures[23-27].

Eqs. (2)–(4) define an extended version of the FK equation, enriched with nonlinear glioma cell diffusion and proliferation terms. The latter terms depend on the oxygen concentration in the tumor microenvironment, which is in turn coupled to cell density through the oxygen consumption term. The functional vascular density controls the supply of oxygen to the tissue. Blood vessel density increases due to tumor angiogenesis and decreases because of vaso-occlusion by high tumor cell density. The values of the parameters used in the simulations and their descriptions are given in Table S1. In addition, a typical full model simulation is shown in Fig. S3 for a representative patient.

We solve the system in Eqs. (2)–(4) by imposing the initial conditions:where the positive parameters c0, n0, and v0 are the initial density of glioma cells spatially distributed in a segment of length ε, the density of functional tumor vasculature, and the oxygen concentration, respectively. Then, L > 0 is the length of the one-dimensional computational domain. In addition, we consider an isolated host tissue in which all system behaviors arise solely due to the interaction terms in Eqs. (2)–(4). This assumption results in no-flux boundary conditions of the form:Both the full and FK models are used to calculate two clinical outputs, namely the tumor IW and TS. The IW at a specific time is defined by the difference between the points where glioma cell density is 80% and 2% of the maximum cellular density. In turn, the TS is obtained by integrating the spatial profile of tumor density and dividing it for the maximum value of the latter.

We run the full model and simulate the growth of the tumor for N patients, each one from a parameter set taken randomly from a uniform distribution over the parameter range. We run simulations for N = 50,100, 250, and 500 with 10 repetitions within each N-case. To generate the patients, we vary five parameters in the list in Table S1, namely the tumor motility D, proliferation rate b, oxygen consumption h1, vascular formation, and occlusion rates g1 and g2, respectively. Then, we use the tumor density at the time of clinical presentation, t0, as the initial condition for the FK model. The latter model is employed to generate predictions at the prediction time t. We also consider the unmodelable variables and clinical outputs at the diagnosis time t, taken randomly between t0 ± 6 months, to build the patient ensemble. Finally, we use the results of the full model in terms of clinical outputs as the ground truth to be compared with the predictions of the FK model alone and with the ones obtained by the BaM3 method.

Probability distribution from the FK model

As described in the previous sections, we take the spatial profile of tumor density at the clinical presentation time t0 as the initial condition of the FK model. We use the latter mathematical model to run simulations over the whole parameter set for cell motility D and proliferation rate b. Then, we define the model-derived pdf as in the following. For each couple of clinical outputs IW* and TS* we calculate the area A(IW*,TS*) over the (IW,TS) plane as A = [(1 − α)IW* < IW < (1 + α)IW*), (1 − α)TS* < TS < (1 + α)TS*)], where α is a given tolerance (here set α = 0.05). Then, we calculate the pdf by normalizing A by the total area of predicted IW and TS values. We store the value of the probability for each patient at the different prediction times and use it to compute the expected value of the model pdf.

Probability distribution of the unmodelables from the full model

To retrieve the data-derived pdf we use a normal kernel density estimator (KDE)[28,29], which depends upon all the data points in the patient ensemble. Briefly, the method estimates the joint probability from which the ensemble entries are drawn through the sum of a kernel function over all the occurrences of the dataset. The kernel function is characterized by a hyperparameter, the bandwidth , which we assume according to Silverman’s rule of thumbwhere d is the number of dimensions, n is the number of observations, and σ is the standard deviation of the ith variate[30]. After calculating , we specify the realization of a specific patient and calculate the value of over the (IW,TS) space of the estimated clinical outputs.

Scoring glioma growth predictions

We calculate for each patient the relative errors d and d as described the main text. To assess how the BaM3 method has changed the prediction of the mathematical model, we compare the latter quantities: if ∣d − d∣ ≤ εd, then there was no change; if d > (1 + ε)d, then the method deteriorated the prediction of the model; if d < (1 − ε)d, then the method improved the prediction of the model. Here, ε is a tolerance used for the comparison, taken to be ε = 0.05.

Calculation of the effective variance

To calculate the effective variance s, we first calculate the mixed central moments Σ of the pdf of interest according to the formulawhere y1 and y2 are the clinical outputs (IW and TS, respectively) and μ1, μ2 the expected values of the corresponding variables. The elements of Σ form a symmetric two-dimensional matrix, for which we calculate the determinant. We define the effective variance s as the natural logarithm of the latter determinant. In Eq. (14) we consider f(y1, y2) to be the pdf from the mathematical model or from the BaM3 method depending on whether we are interested in the effective variance s or s, respectively.

Pdf from the two-compartment model in CLL

Messmer and colleagues[31] measured the fraction of labeled B-CLL cells in a cohort of 17 CLL patients that were administered deuterated water. They calibrated a two-compartment model on each patient and were able to reproduce the kinetics of labeled cells over a long time. We adopt their model and use it to generate the pdf for the CLL example. The fraction of labeled cells over time is calculated through the expressionwhere g(0) is the initial fraction of cells in the first compartment, b the fractional cell birth, v the relative size of the compartments, and h(t) is the deuterated water concentration of the body over time. The latter is a function of the fractional daily water exchange f. We refer the interested reader to the supplementary information of Messmer et al.[31] for a more detailed description of the model and a full account of the model parameters. In this work we focus on three quantities, namely b, v, and f, and run the model in Eq. (15) over the experimental range. This range was obtained by considering the patient-specific fitting performed by Messmer and colleagues and selecting the minimum and maximum values. We evaluate the fraction of labeled cells at day 50, f50, and build the probability distribution from its histogram, by counting the number of occurrences of a given for and then normalizing the result. For the CLL example, all the patients start with the same initial fraction of labeled cells, set to zero.

Pdf from the patients’ unmodelables in CLL

The data-derived pdf in the CLL example is obtained from four unmodelable quantities that are measured for each patient during the study. We consider all the possible combinations of unmodelables and calculate the mean squared error (MSE) for each case. The scatter plot in the same picture refers to the case in which the CD38 expression (x), age (x), growth rate of white blood cells (x), and VH mutation status (x) are added consecutively with the specified order. As in the glioma example, we build the sub-dataset (y, x), where y and x = (x) are the f50 and unmodelable variables of each patient, respectively, and apply the KDE using Silverman’s rule for the hyperparameters. The requested pdf, i.e., is obtained by conditioning the probability from the KDE with the realizations of the unmodelables of the specific patient and calculating the result over the range of the estimated clinical output .

Mathematical model for ovarian cancer

We assume the total number of tumor cells T to be composed of the sensitive S and resistant R subpopulations. The latter are described by the following system of ordinary differential equations (ODEs):where γ is the tumor net growth rate, δ = δ(t) is the death rate induced by chemotherapy, τ is the mutation rate from sensitive to resistant cells, and λ is a factor that accounts for reduced death by therapy in resistant cells. As detailed in Fig. S12, the treatment is composed of three phases: first, the patients undergo different cycles of NACT; then, surgery is performed. The latter reduces the total tumor volume, irrespective of cells being sensitive or resistant, of a factor β. Finally, another series of chemotherapy cycles is performed. During chemotherapy, δ = δ0, whereas we set this parameter to zero after chemotherapy and until tumor relapse. The latter condition occurs when T reaches the value T.

Equations (16) and (17) can be analytically integrated, and their results used to build the probability distribution of the clinical output—TtR, in this case. To obtain the pdf from the model, we calculate the time the tumor takes to reach the cell number at relapse T starting from the cell number after therapy. We perform this calculation using the initial tumor cell number of each patient, and by varying both the initial fraction of S cells, x0, and the chemotherapy-induced death rate, δ0. We then obtain the patient-specific probability distribution from the histogram of TtR, similarly to what is done in the previous section for CLL. For x0, we select a range between 0.4 and 0.9, accounting for tumors with different initial degrees of intrinsic resistance[32]. For δ0, we first use a uniform distribution between 0.1 and 10 days−1, accounting for a wide variation in death rates. The latter choice produces an almost flat distribution for the clinical output (see Fig. S13). To improve the mathematical model parametrization, we use the information about the tumor volume change after the first cycle of chemotherapy, which is included in the dataset. By fitting T obtained from Eqs. (16) and (17) to the observed volume change, we find a value of δ0 for each patient in the dataset[32]. We take the mean value of these rates and use it to update the model pdf (see Fig. S14). We consider a range for δ0 that is centered around its mean value across the patients, within an interval of ±40%. Selecting other ranges provides similar results, however, a variation of 40% returns the lowest MSE. Analytical integration of Eqs. (16) and (17), as well as additional details about model parametrization are available in Supplementary Note 2.

Unmodelable variable for the ovarian cancer study

We build the data-derived pdf for the ovarian cancer example by exploiting the information about the age of the patients at diagnosis. Similarly to what done in the previous test cases, we first build the sub-dataset (TtR, A) by entering the information of each patient (here, A is the patient age). Then, we apply the KDE using Silverman’s rule to estimate the bandwidth and calculate the joint probability . The data-derived pdf for each patient is finally obtained over the domain of the clinical output TtR by considering the patient-specific age A = A*.