Global stability and parameter analysis reinforce therapeutic targets of PD-L1-PD-1 and MDSCs for glioblastoma.

Glioblastoma (GBM) is an aggressive primary brain cancer that currently has minimally effective treatments. Like other cancers, immunosuppression by the PD-L1-PD-1 immune checkpoint complex is a prominent axis by which glioma cells evade the immune system. Myeloid-derived suppressor cells (MDSCs), which are recruited to the glioma microenviroment, also contribute to the immunosuppressed GBM microenvironment by suppressing T cell functions. In this paper, we propose a GBM-specific tumor-immune ordinary differential equations model of glioma cells, T cells, and MDSCs to provide theoretical insights into the interactions between these cells. Equilibrium and stability analysis indicates that there are unique tumorous and tumor-free equilibria which are locally stable under certain conditions. Further, the tumor-free equilibrium is globally stable when T cell activation and the tumor kill rate by T cells overcome tumor growth, T cell inhibition by PD-L1-PD-1 and MDSCs, and the T cell death rate. Bifurcation analysis suggests that a treatment plan that includes surgical resection and therapeutics targeting immune suppression caused by the PD-L1-PD1 complex and MDSCs results in the system tending to the tumor-free equilibrium. Using a set of preclinical experimental data, we implement the approximate Bayesian computation (ABC) rejection method to construct probability density distributions that estimate model parameters. These distributions inform an appropriate search curve for global sensitivity analysis using the extended fourier amplitude sensitivity test. Sensitivity results combined with the ABC method suggest that parameter interaction is occurring between the drivers of tumor burden, which are the tumor growth rate and carrying capacity as well as the tumor kill rate by T cells, and the two modeled forms of immunosuppression, PD-L1-PD-1 immune checkpoint and MDSC suppression of T cells. Thus, treatment with an immune checkpoint inhibitor in combination with a therapeutic targeting the inhibitory mechanisms of MDSCs should be explored.

Dynamics of a linearly perturbed May-Leonard competition model.

The May-Leonard model was introduced to examine the behavior of three competing populations where rich dynamics, such as limit cycles and nonperiodic cyclic solutions, arise. In this work, we perturb the system by adding the capability of global mutations, allowing one species to evolve to the other two in a linear manner. We find that for small mutation rates, the perturbed system not only retains some of the dynamics seen in the classical model, such as the three-species equal-population equilibrium bifurcating to a limit cycle, but also exhibits new behavior. For instance, we capture curves of fold bifurcations where pairs of equilibria emerge and then coalesce. As a result, we uncover parameter regimes with new types of stable fixed points that are distinct from the single- and dual-population equilibria characteristic of the original model. On the contrary, the linearly perturbed system fails to maintain heteroclinic connections that exist in the original system. In short, a linear perturbation proves to be significant enough to substantially influence the dynamics, even with small mutation rates.

An Agent-Based Model of Biting Midge Dynamics to Understand Bluetongue Outbreaks.

Bluetongue (BT) is a well-known vector-borne disease that infects ruminants such as sheep, cattle, and deer with high mortality rates. Recent outbreaks in Europe highlight the importance of understanding vector-host dynamics and potential courses of action to mitigate the damage that can be done by BT. We present an agent-based model, entitled 'MidgePy', that focuses on the movement of individual Culicoides spp. biting midges and their interactions with ruminants to understand their role as vectors in BT outbreaks, especially in regions that do not regularly experience outbreaks. The results of our sensitivity analysis suggest that midge survival rate has a significant impact on the probability of a BTV outbreak as well as its severity. Using midge flight activity as a proxy for temperature, we found that an increase in environmental temperature corresponded with an increased probability of outbreak after identifying parameter regions where outbreaks are more likely to occur. This suggests that future methods to control BT spread could combine large-scale vaccination programs with biting midge population control measures such as the use of pesticides. Spatial heterogeneity in the environment is also explored to give insight on optimal farm layouts to reduce the potential for BT outbreaks.

Global stability and parameter analysis reinforce therapeutic targets of PD-L1-PD-1 and MDSCs for glioblastoma.

Glioblastoma (GBM) is an aggressive primary brain cancer that currently has minimally effective treatments. Like other cancers, immunosuppression by the PD-L1-PD-1 immune checkpoint complex is a prominent axis by which glioma cells evade the immune system. Myeloid-derived suppressor cells (MDSCs), which are recruited to the glioma microenviroment, also contribute to the immunosuppressed GBM microenvironment by suppressing T cell functions. In this paper, we propose a GBM-specific tumor-immune ordinary differential equations model of glioma cells, T cells, and MDSCs to provide theoretical insights into the interactions between these cells. Equilibrium and stability analysis indicates that there are unique tumorous and tumor-free equilibria which are locally stable under certain conditions. Further, the tumor-free equilibrium is globally stable when T cell activation and the tumor kill rate by T cells overcome tumor growth, T cell inhibition by PD-L1-PD-1 and MDSCs, and the T cell death rate. Bifurcation analysis suggests that a treatment plan that includes surgical resection and therapeutics targeting immune suppression caused by the PD-L1-PD1 complex and MDSCs results in the system tending to the tumor-free equilibrium. Using a set of preclinical experimental data, we implement the Approximate Bayesian Computation (ABC) rejection method to construct probability density distributions that estimate model parameters. These distributions inform an appropriate search curve for global sensitivity analysis using the extended Fourier Amplitude Sensitivity Test (eFAST). Sensitivity results combined with the ABC method suggest that parameter interaction is occurring between the drivers of tumor burden, which are the tumor growth rate and carrying capacity as well as the tumor kill rate by T cells, and the two modeled forms of immunosuppression, PD-L1-PD-1 immune checkpoint and MDSC suppression of T cells. Thus, treatment with an immune checkpoint inhibitor in combination with a therapeutic targeting the inhibitory mechanisms of MDSCs should be explored.

Spreading mechanics and differentiation of astrocytes during retinal development.

The retinal vasculature supplies oxygen to the inner layers of the retina, the light-sensitive tissue in the eye. During development, formation of the retinal vasculature depends on prior establishment of a mesh of astrocytes, a type of glial cell, which guide the growth of the vascular network. Astrocytes emerge from the optic nerve head and proliferate and spread, forming a mesh-like layer over the retinal surface. The initially formed cells are termed astrocyte precursor cells (APCs), which differentiate into immature perinatal astrocytes (IPAs) during the prenatal period. A continuum model is developed to describe the proliferation, differentiation, and migration these cells. Effects of oxygen and growth factor levels on proliferation and differentiation are included. Cell migration is driven by gradients in tension in the astrocyte mesh, which varies inversely with total density. The resulting governing equations have the form of a nonlinear diffusion-like equation. The model can account for the observed radial spread over time of the astrocyte disk. Experimental observations show that the APCs form a narrow rim around the edge of this disk, with IPAs in the interior. The model predicts this behavior if the mobility of the APCs is assumed to be higher than that of the IPAs under a given tension gradient. Thus, the model shows how tension-driven cell motions can account for separation of cell types in a cell layer spreading over a substrate.

Using a continuum model to decipher the mechanics of embryonic tissue spreading from time-lapse image sequences: An approximate Bayesian computation approach.

Advanced imaging techniques generate large datasets capable of describing the structure and kinematics of tissue spreading in embryonic development, wound healing, and the progression of many diseases. These datasets can be integrated with mathematical models to infer biomechanical properties of the system, typically identifying an optimal set of parameters for an individual experiment. However, these methods offer little information on the robustness of the fit and are generally ill-suited for statistical tests of multiple experiments. To overcome this limitation and enable efficient use of large datasets in a rigorous experimental design, we use the approximate Bayesian computation rejection algorithm to construct probability density distributions that estimate model parameters for a defined theoretical model and set of experimental data. Here, we demonstrate this method with a 2D Eulerian continuum mechanical model of spreading embryonic tissue. The model is tightly integrated with quantitative image analysis of different sized embryonic tissue explants spreading on extracellular matrix (ECM) and is regulated by a small set of parameters including forces on the free edge, tissue stiffness, strength of cell-ECM adhesions, and active cell shape changes. We find statistically significant trends in key parameters that vary with initial size of the explant, e.g., for larger explants cell-ECM adhesion forces are weaker and free edge forces are stronger. Furthermore, we demonstrate that estimated parameters for one explant can be used to predict the behavior of other similarly sized explants. These predictive methods can be used to guide further experiments to better understand how collective cell migration is regulated during development.

Mathematical Analysis of Glioma Growth in a Murine Model.

Five immunocompetent C57BL/6-cBrd/cBrd/Cr (albino C57BL/6) mice were injected with GL261-luc2 cells, a cell line sharing characteristics of human glioblastoma multiforme (GBM). The mice were imaged using magnetic resonance (MR) at five separate time points to characterize growth and development of the tumor. After 25 days, the final tumor volumes of the mice varied from 12 mm3 to 62 mm3, even though mice were inoculated from the same tumor cell line under carefully controlled conditions. We generated hypotheses to explore large variances in final tumor size and tested them with our simple reaction-diffusion model in both a 3-dimensional (3D) finite difference method and a 2-dimensional (2D) level set method. The parameters obtained from a best-fit procedure, designed to yield simulated tumors as close as possible to the observed ones, vary by an order of magnitude between the three mice analyzed in detail. These differences may reflect morphological and biological variability in tumor growth, as well as errors in the mathematical model, perhaps from an oversimplification of the tumor dynamics or nonidentifiability of parameters. Our results generate parameters that match other experimental in vitro and in vivo measurements. Additionally, we calculate wave speed, which matches with other rat and human measurements.

A data-motivated density-dependent diffusion model of in vitro glioblastoma growth.

Glioblastoma multiforme is an aggressive brain cancer that is extremely fatal. It is characterized by both proliferation and large amounts of migration, which contributes to the difficulty of treatment. Previous models of this type of cancer growth often include two separate equations to model proliferation or migration. We propose a single equation which uses density-dependent diffusion to capture the behavior of both proliferation and migration. We analyze the model to determine the existence of traveling wave solutions. To prove the viability of the density-dependent diffusion function chosen, we compare our model with well-known in vitro experimental data.

TGF-mediated dynamics in a system of many coupled nephrons.

This paper presents a mathematical model of a system of many coupled nephrons branching from a common cortical radial artery, and accompanying analysis of that system. This modeling effort is a first step in understanding how coupling magnifies the tendency of nephrons to oscillate owing to tubuloglomerular feedback. Central to the present work is the single nephron integral model (as in Pitman et al., The IMA Volumes in Mathematics and Its Applications, vol. 129, pp. 345-364, 2002 and in Zaritski, Ph.D. Dissertation, 1999) which is a simplification of the single nephron PDE model of Layton et al. (Am. J. Physiol. 261, F904-F919, 1991). A second principal idea used in the present model is a coupling of model nephrons, generalizing the work of Pitman et al. (Bull. Math. Biol. 66, 1463-1492, 2004) who proposed a model of two coupled nephrons. In this study, we couple nephrons through a nearest neighbor interaction.Speaking generally, our results suggest that a series of similar nephrons coupled to their nearest neighbors are more prone to be found in an oscillatory mode, relative to a single nephron with the same properties. More specifically, we show analytically that, for N coupled identical nephrons, the region supporting oscillatory solutions in the time delay-gain parameter plane increases with N. Numerical simulations suggest that, if N nephrons have gains and time delays that do not differ by much, the system is, again, more prone to oscillate, relative to a single nephron, and the oscillations tend to be approximately synchronous and in-phase. We examine the effect of parameters on bifurcation. We also examine alternative models of coupling; this analysis allows us to conclude that the increased propensity of coupled nephrons to oscillate is a robust finding, true for several models of nephron interaction.