Classical mathematical models for prediction of response to chemotherapy and immunotherapy

Classical mathematical models of tumor growth have shaped our understanding of cancer and have broad practical implications for treatment scheduling and dosage. However, even the simplest textbook models have been barely validated in real world-data of human patients. In this study, we fitted a range of differential equation models to tumor volume measurements of patients undergoing chemotherapy or cancer immunotherapy for solid tumors. We used a large dataset of 1472 patients with three or more measurements per target lesion, of which 652 patients had six or more data points. We show that the early treatment response shows only moderate correlation with the final treatment response, demonstrating the need for nuanced models. We then perform a head-to-head comparison of six classical models which are widely used in the field: the Exponential, Logistic, Classic Bertalanffy, General Bertalanffy, Classic Gompertz and General Gompertz model. Several models provide a good fit to tumor volume measurements, with the Gompertz model providing the best balance between goodness of fit and number of parameters. Similarly, when fitting to early treatment data, the general Bertalanffy and Gompertz models yield the lowest mean absolute error to forecasted data, indicating that these models could potentially be effective at predicting treatment outcome. In summary, we provide a quantitative benchmark for classical textbook models and state-of-the art models of human tumor growth. We publicly release an anonymized version of our original data, providing the first benchmark set of human tumor growth data for evaluation of mathematical models.     

A comparison of three models of the inositol trisphosphate receptor

The inositol (1,4,5)-trisphosphate receptor (IPR) plays a crucial role in calcium dynamics in a wide range of cell types, and is often a central feature in quantitative models of calcium oscillations and waves. We compare three mathematical models of the IPR, fitting each of them to the same data set to determine ranges for the parameter values. Each of the fits indicates that fast activation of the receptor, followed by slow inactivation, is an important feature of the model, and also that the speed of inositol trisphosphate (IP₃) binding cannot necessarily be assumed to be faster than Ca²⁺ activation. In addition, the model which assumed saturating binding rates of Ca²⁺ to the IPR demonstrated the best fit. However, lack of convergence in the fitting procedure indicates that responses to step increases of [Ca²⁺] and [IP₃] provide insufficient data to determine the parameters unambiguously in any of the models.     

When the optimal is not the best: parameter estimation in complex biological models

Background

The vast computational resources that became available during the past decade enabled the development and simulation of increasingly complex mathematical models of cancer growth. These models typically involve many free parameters whose determination is a substantial obstacle to model development. Direct measurement of biochemical parameters in vivo is often difficult and sometimes impracticable, while fitting them under data-poor conditions may result in biologically implausible values.

Results

We discuss different methodological approaches to estimate parameters in complex biological models. We make use of the high computational power of the Blue Gene technology to perform an extensive study of the parameter space in a model of avascular tumor growth. We explicitly show that the landscape of the cost function used to optimize the model to the data has a very rugged surface in parameter space. This cost function has many local minima with unrealistic solutions, including the global minimum corresponding to the best fit.

Conclusions

The case studied in this paper shows one example in which model parameters that optimally fit the data are not necessarily the best ones from a biological point of view. To avoid force-fitting a model to a dataset, we propose that the best model parameters should be found by choosing, among suboptimal parameters, those that match criteria other than the ones used to fit the model. We also conclude that the model, data and optimization approach form a new complex system and point to the need of a theory that addresses this problem more generally.     

Markov chain Monte Carlo fitting of single-channel data from inositol trisphosphate receptors

In many cell types, the inositol trisphosphate receptor (IPR) is one of the important components that control intracellular calcium dynamics, and an understanding of this receptor (which is also a calcium channel) is necessary for an understanding of calcium oscillations and waves. Recent advances in experimental techniques now allow for the measurement of single-channel activity of the IPR in conditions similar to its native environment, and these data can be used to determine the rate constants in Markov models of the IPR. We illustrate a parameter estimation method based on Markov chain Monte Carlo, which can be used to fit directly to single-channel data, and determining, as an intrinsic part of the fit, the times at which the IPR is opening and closing. We show, using simulated data, the most complex Markov model that can be unambiguously determined from steady-state data and show that non-steady-state data is required to determine more complex models.     

Positive feedback and angiogenesis in tumor growth control

In vivo tumor growth data from experiments performed in our laboratory suggest that basic fibroblast growth factor (bFGF) and vascular endothelial growth factor (VEGF) are angiogenic signals emerging from an up-regulated genetic message in the proliferating rim of a solid tumor in response to tumor-wide hypoxia. If these signals are generated in response to unfavorable environmental conditions, i.e. a decrease in oxygen tension, then the tumor may play an active role in manipulating its own environment. We have idealized this type of adaptive behavior in our mathematical model via a parameter which represents the carrying capacity of the host for the tumor. If that model parameter is held constant, then environmental control is limited to tumor shape and mitogenic signal processing. However, if we assume that the response of the local stroma to these signals is an increase in the host's ability to support an ever larger tumor, then our models describe a positive feedback control system. In this paper, we generalize our previous results to a model including a carrying capacity which depends on the size of the proliferating compartment in the tumor. Specific functional forms for the carrying capacity are discussed. Stability criteria of the system and steady state conditions for these candidate functions are analyzed. The dynamics needed to generate stable tumor growth, including countervailing negative feedback signals, are discussed in detail with respect to both their mathematical and biological properties.     

Mathematical and numerical comparisons of five single-population growth models

We investigate the properties of five mathematical models used to represent the growth of a single population. By imposing a common set of (normalizing) initial conditions, we are able to calculate and explicitly compare the time intervals required to reach specific values of population levels. Based on these results, we conclude that one must be careful when applying these models to interpret the dynamics of single-population growth. An additional implication is that they provide evidence that such caution should also be extended to the incorporation of these models into the formulation of interacting, multi-population models, which are used, for example, to study the spread of disease.     

Comparing the consequences of induced and constitutive plant resistance for herbivore population dynamics

Although it has been suggested that induced and constitutive plant resistance should have different effects on insect herbivore population dynamics, there is little experimental evidence that plant resistance can influence herbivore populations longer than one season. We used a density-manipulation experiment and model fitting to examine the effects of constitutive and induced resistance on herbivore dynamics over both the short and long term. We used likelihood methods to fit population dynamic models to recruitment data for populations of Mexican bean beetles on soybean varieties with no resistance, constitutive resistance, or induced resistance. We compared model configurations that fit parameters for resistance types separately to models that did not account for resistance type. Models representing the hypothesis that the three resistance types differed in their effects on beetle dynamics received the most support. Induced resistance resulted in lower population growth rates and stronger density dependence than no resistance. Constitutive resistance resulted in lower population growth rates and stronger density dependence than induced resistance. Constitutive resistance had a stronger effect on both short-term beetle recruitment and predicted beetle population dynamics than induced resistance. The results of this study suggest that induced and constitutive resistance can differ in their effects on herbivore populations even in a relatively complex system.     

Asymmetric regression models with limited responses with an application to antibody response to vaccine

We develop regression models for limited and censored data based on the mixture between the log‐power‐normal and Bernoulli‐type distributions. A likelihood‐based approach is implemented for parameter estimation and a small‐scale simulation study is conducted to evaluate parameter recovery, with emphasis on bias estimation. The main conclusion is that the approach is very much satisfactory for moderate and large sample sizes. A real data example, the safety and immunogenecity study of measles vaccine in Haiti, is presented to illustrate how different models can be used to fit this type of data. As shown, the asymmetric models considered seem to present the best fit for the data set under study, revealing significance of the explanatory variable sex, which is not found significant with the log‐normal model.     

Keratin Dynamics: Modeling the Interplay between Turnover and Transport

Keratin are among the most abundant proteins in epithelial cells. Functions of the keratin network in cells are shaped by their dynamical organization. Using a collection of experimentally-driven mathematical models, different hypotheses for the turnover and transport of the keratin material in epithelial cells are tested. The interplay between turnover and transport and their effects on the keratin organization in cells are hence investigated by combining mathematical modeling and experimental data. Amongst the collection of mathematical models considered, a best model strongly supported by experimental data is identified. Fundamental to this approach is the fact that optimal parameter values associated with the best fit for each model are established. The best candidate among the best fits is characterized by the disassembly of the assembled keratin material in the perinuclear region and an active transport of the assembled keratin. Our study shows that an active transport of the assembled keratin is required to explain the experimentally observed keratin organization.     

Development of mathematical models (Logistic, Gompertz and Richards models) describing the growth pattern of Pseudomonas putida (NICM 2174)

Bacterial growth curve, which is asymptotic after a certain period, is described using three different mathematical models, namely, Logistic model, Gompertz model and Richards model. The equations for these three models are fitted by evaluating the mathematical parameters involved in these models. This is done by applying the method of partial sums to the data in Table 1 containing the optical density values for different cell mass at different time intervals. The sum of square of residuals between the expected optical density values and the experimental values is calculated for each of these models. In the cases tested, the Logistic model was found to be the best fit for the growth curve of Pseudomonas putida (NICM 2174) and was found to be easy to use. These results fit the data very well at 5% level for more than 70% of the readings.     

A model of HIV-1 infection with two time delays: mathematical analysis and comparison with patient data

Mathematical models have made considerable contributions to our understanding of HIV dynamics. Introducing time delays to HIV models usually brings challenges to both mathematical analysis of the models and comparison of model predictions with patient data. In this paper, we incorporate two delays, one the time needed for infected cells to produce virions after viral entry and the other the time needed for the adaptive immune response to emerge to control viral replication, into an HIV-1 model. We begin model analysis with proving the positivity and boundedness of the solutions, local stability of the infection-free and infected steady states, and uniform persistence of the system. By developing a few Lyapunov functionals, we obtain conditions ensuring global stability of the steady states. We also fit the model including two delays to viral load data from 10 patients during primary HIV-1 infection and estimate parameter values. Although the delay model provides better fits to patient data (achieving a smaller error between data and modeling prediction) than the one without delays, we could not determine which one is better from the statistical standpoint. This highlights the need of more data sets for model verification and selection when we incorporate time delays into mathematical models to study virus dynamics.     

A model for hormonal control of the menstrual cycle: Structural consistency but sensitivity with regard to data

This study presents a 13-dimensional system of delayed differential equations which predicts serum concentrations of five hormones important for regulation of the menstrual cycle. Parameters for the system are fit to two different data sets for normally cycling women. For these best fit parameter sets, model simulations agree well with the two different data sets but one model also has an abnormal stable periodic solution, which may represent polycystic ovarian syndrome. This abnormal cycle occurs for the model in which the normal cycle has estradiol levels at the high end of the normal range. Differences in model behavior are explained by studying hysteresis curves in bifurcation diagrams with respect to sensitive model parameters. For instance, one sensitive parameter is indicative of the estradiol concentration that promotes pituitary synthesis of a large amount of luteinizing hormone, which is required for ovulation. Also, it is observed that models with greater early follicular growth rates may have a greater risk of cycling abnormally.     

Density-structured models for plant population dynamics

Density-structured models are structured population models in which the state variable is the proportion of populations or sites in a small number of discrete density states. Although such models have rarely been used, they have the advantage that they are straightforward to parameterize, make few assumptions about population dynamics, and permit rapid data collection using coarse density assessment. In this article, we highlight their use in relating population dynamics to environmental variation and their robustness to measurement error. We show that density-structured models are able to accurately represent population dynamics under a wide range of conditions. We look at the effects of including a persistent seedbank and describe numerical approximations for the mean and variance of population size. For simulated data, we determine the extent to which the underlying continuous process may be inferred from density-structured data. Finally, we discuss issues of parameter estimation and applications for which these types of models may be useful.     

Towards Predictive Computational Models of Oncolytic Virus Therapy: Basis for Experimental Validation and Model Selection

Oncolytic viruses are viruses that specifically infect cancer cells and kill them, while leaving healthy cells largely intact. Their ability to spread through the tumor makes them an attractive therapy approach. While promising results have been observed in clinical trials, solid success remains elusive since we lack understanding of the basic principles that govern the dynamical interactions between the virus and the cancer. In this respect, computational models can help experimental research at optimizing treatment regimes. Although preliminary mathematical work has been performed, this suffers from the fact that individual models are largely arbitrary and based on biologically uncertain assumptions. Here, we present a general framework to study the dynamics of oncolytic viruses that is independent of uncertain and arbitrary mathematical formulations. We find two categories of dynamics, depending on the assumptions about spatial constraints that govern that spread of the virus from cell to cell. If infected cells are mixed among uninfected cells, there exists a viral replication rate threshold beyond which tumor control is the only outcome. On the other hand, if infected cells are clustered together (e.g. in a solid tumor), then we observe more complicated dynamics in which the outcome of therapy might go either way, depending on the initial number of cells and viruses. We fit our models to previously published experimental data and discuss aspects of model validation, selection, and experimental design. This framework can be used as a basis for model selection and validation in the context of future, more detailed experimental studies. It can further serve as the basis for future, more complex models that take into account other clinically relevant factors such as immune responses.     

ODE models for oncolytic virus dynamics

Replicating oncolytic viruses are able to infect and lyse cancer cells and spread through the tumor, while leaving normal cells largely unharmed. This makes them potentially useful in cancer therapy, and a variety of viruses have shown promising results in clinical trials. Nevertheless, consistent success remains elusive and the correlates of success have been the subject of investigation, both from an experimental and a mathematical point of view. Mathematical modeling of oncolytic virus therapy is often limited by the fact that the predicted dynamics depend strongly on particular mathematical terms in the model, the nature of which remains uncertain. We aim to address this issue in the context of ODE modeling, by formulating a general computational framework that is independent of particular mathematical expressions. By analyzing this framework, we find some new insights into the conditions for successful virus therapy. We find that depending on our assumptions about the virus spread, there can be two distinct types of dynamics. In models of the first type (the “fast spread” models), we predict that the viruses can eliminate the tumor if the viral replication rate is sufficiently high. The second type of models is characterized by a suboptimal spread (the “slow spread” models). For such models, the simulated treatment may fail, even for very high viral replication rates. Our methodology can be used to study the dynamics of many biological systems, and thus has implications beyond the study of virus therapy of cancers.     

Action of S-Carbamoyl and S-Thiocarbamoyl Derivatives of L-Cysteine on L-Amino Acid Oxidase

Methods are investigated for evaluating the kinetic parameters in a modified Monod’s equation which give the best fit to the growth thermograms for bacterial cultures observed in batch calorimeters. Four mathematical methods were employed as parameter fitting techniques. The growth thermograms observed for soil microbes cultured with glucose as a limiting substrate were used as the objects of the analysis. For the calculation of the heat evolution rate, the Runge-Kutta method, which is commonly used for the numerical analysis, was employed. A comparison of the results obtained by the four methods in terms of closeness of fit to the actual thermograms showed that optimization by direct searching with the Simplex method is the most effective procedure for obtaining the best values of the parameters to reproduce the observed thermograms.     

A glucose-insulin pharmacodynamic surface modeling validation and comparison of metabolic system models

Metabolic system modeling for model-based glycaemic control is becoming increasingly important. Few metabolic system models are clinically validated for both fit to the data and prediction ability. This research introduces a new additional form of pharmaco-dynamic (PD) surface comparison for model analysis and validation. These 3D surfaces are developed for 3 clinically validated models and 1 model with an added saturation dynamic. The models include the well-known Minimal Model. They are fit to two different data sets of clinical PD data from hyperinsulinaemic clamp studies at euglycaemia and/or hyperglycaemia. The models are fit to the first data set to determine an optimal set of population parameters. The second data set is used to test trend prediction of the surface modeling as it represents a lower insulin sensitivity cohort and should thus require only scaling in these (or related) parameters to match this data set. This particular approach clearly highlights differences in modeling methods, and the model dynamics utilized that may not appear as clearly in other fitting or prediction validation methods.Across all models saturation of insulin action is seen to be an important determinant of prediction and fit quality. In particular, the well-reported under-modeling of insulin sensitivity in the Minimal Model can be seen in this context to be a result of a lack of saturation dynamics, which in turn affects its ability to detect differences between cohorts. The overall approach of examining PD surfaces is seen to be an effective means of analyzing and thus validating a metabolic model's inherent dynamics and basic trend prediction on a population level, but is not a replacement for data driven, patient-specific fit and prediction validation for clinical use. The overall method presented could be readily generalized to similar PD systems and therapeutics.     

Comparison of amplitude recovery dynamics of two limit cycle oscillator models of the human circadian pacemaker

At an organism level, the mammalian circadian pacemaker is a two-dimensional system. For these two dimensions, phase (relative timing) and amplitude of the circadian pacemaker are commonly used. Both the phase and the amplitude (A) of the human circadian pacemaker can be observed within multiple physiological measures--including plasma cortisol, plasma melatonin, and core body temperature (CBT)--all of which are also used as markers of the circadian system. Although most previous work has concentrated on changes in phase of the circadian system, critically timed light exposure can significantly reduce the amplitude of the pacemaker. The rate at which the amplitude recovers to its equilibrium level after reduction can have physiological significance. Two mathematical models that describe the phase and amplitude dynamics of the pacemaker have been reported. These models are essentially equivalent in predictions of phase and in predictions of amplitude recovery for small changes from an equilibrium value (A = 1), but are markedly different in the prediction of recovery rates when A < 0.6. To determine which dynamic model best describes the amplitude recovery observed in experimental data; both models were fit to CBT data using a maximum likelihood procedure and compared using Akaike's Information Criterion (AIC). For all subjects, the model with the lower recovery rate provided a better fit to data in terms of AIC, supporting evidence that the amplitude recovery of the endogenous pacemaker is slow at low amplitudes. Experiments derived from model predictions are proposed to test the influence of low amplitude recovery on the physiological and neurobehavioral functions.