A Mathematical Model of TCR-T Cell Therapy for Cervical Cancer.

Engineered T cell receptor (TCR)-expressing T (TCR-T) cells are intended to drive strong anti-tumor responses upon recognition of the specific cancer antigen, resulting in rapid expansion in the number of TCR-T cells and enhanced cytotoxic functions, causing cancer cell death. However, although TCR-T cell therapy against cancers has shown promising results, it remains difficult to predict which patients will benefit from such therapy. We develop a mathematical model to identify mechanisms associated with an insufficient response in a mouse cancer model. We consider a dynamical system that follows the population of cancer cells, effector TCR-T cells, regulatory T cells (Tregs), and "non-cancer-killing" TCR-T cells. We demonstrate that the majority of TCR-T cells within the tumor are "non-cancer-killing" TCR-T cells, such as exhausted cells, which contribute little or no direct cytotoxicity in the tumor microenvironment (TME). We also establish two important factors influencing tumor regression: the reversal of the immunosuppressive TME following depletion of Tregs, and the increased number of effector TCR-T cells with antitumor activity. Using mathematical modeling, we show that certain parameters, such as increasing the cytotoxicity of effector TCR-T cells and modifying the number of TCR-T cells, play important roles in determining outcomes.

Minimal Mechanisms of Microtubule Length Regulation in Living Cells.

The microtubule cytoskeleton is responsible for sustained, long-range intracellular transport of mRNAs, proteins, and organelles in neurons. Neuronal microtubules must be stable enough to ensure reliable transport, but they also undergo dynamic instability, as their plus and minus ends continuously switch between growth and shrinking. This process allows for continuous rebuilding of the cytoskeleton and for flexibility in injury settings. Motivated by in vivo experimental data on microtubule behavior in Drosophila neurons, we propose a mathematical model of dendritic microtubule dynamics, with a focus on understanding microtubule length, velocity, and state-duration distributions. We find that limitations on microtubule growth phases are needed for realistic dynamics, but the type of limiting mechanism leads to qualitatively different responses to plausible experimental perturbations. We therefore propose and investigate two minimally-complex length-limiting factors: limitation due to resource (tubulin) constraints and limitation due to catastrophe of large-length microtubules. We combine simulations of a detailed stochastic model with steady-state analysis of a mean-field ordinary differential equations model to map out qualitatively distinct parameter regimes. This provides a basis for predicting changes in microtubule dynamics, tubulin allocation, and the turnover rate of tubulin within microtubules in different experimental environments.

Optimal Rotation Age in Fast Growing Plantations: A Dynamical Optimization Problem.

Forest plantations are economically and environmentally relevant, as they play a key role in timber production and carbon capture. It is expected that the future climate change scenario affects forest growth and modify the rotation age for timber production. However, mathematical models on the effect of climate change on the rotation age for timber production remain still limited. We aim to determine the optimal rotation age that maximizes the net economic benefit of timber volume in a negative scenario from the climatic point of view. For this purpose, a bioeconomic optimal control problem was formulated from a system of Ordinary Differential Equations (ODEs) governed by the state variables live biomass volume, intrinsic growth rate, and area affected by fire. Then, four control variables were associated to the system, representing forest management activities, which are felling, thinning, reforestation, and fire prevention. The existence of optimal control solutions was demonstrated, and the solutions of the optimal control problem were also characterized using Pontryagin's Maximum Principle. The solutions of the model were approximated numerically by the Forward-Backward Sweep method. To validate the model, two scenarios were considered: a realistic scenario that represents current forestry activities for the exotic species Pinus radiata D. Don, and a pessimistic scenario, which considers environmental conditions conducive to a higher occurrence of forest fires. The optimal solution that maximizes the net benefit of timber volume consists of a strategy that considers all four control variables simultaneously. For felling and thinning, regardless of the scenario considered, the optimal strategy is to spend on both activities depending on the amount of biomass in the field. Similarly, for reforestation, the optimal strategy is to spend as the forest is harvested. In the case of fire prevention, in the realistic scenario, the optimal strategy consists of reducing the expenses in fire prevention because the incidence of fires is lower, whereas in the pessimistic scenario, the opposite is true. It is concluded that the optimal rotation age that maximizes the net economic benefit of timber volume in P. radiata plantations is 24 and 19 years for the realistic and pessimistic scenarios, respectively. This corroborates that the presence of fires influences the determination of the optimal rotation age, and as a consequence, the net economic benefit.

Integrating Diversity, Equity, and Inclusion into Preclinical, Clinical, and Public Health Mathematical Models.

Mathematical modelling applied to preclinical, clinical, and public health research is critical for our understanding of a multitude of biological principles. Biology is fundamentally heterogeneous, and mathematical modelling must meet the challenge of variability head on to ensure the principles of diversity, equity, and inclusion (DEI) are integrated into quantitative analyses. Here we provide a follow-up perspective on the DEI plenary session held at the 2023 Society for Mathematical Biology Annual Meeting to discuss key issues for the increased integration of DEI in mathematical modelling in biology.

Dynamics of Information Flow and Task Allocation of Social Insect Colonies: Impacts of Spatial Interactions and Task Switching.

Models of social interaction dynamics have been powerful tools for understanding the efficiency of information spread and the robustness of task allocation in social insect colonies. How workers spatially distribute within the colony, or spatial heterogeneity degree (SHD), plays a vital role in contact dynamics, influencing information spread and task allocation. We used agent-based models to explore factors affecting spatial heterogeneity and information flow, including the number of task groups, variation in spatial arrangements, and levels of task switching, to study: (1) the impact of multiple task groups on SHD, contact dynamics, and information spread, and (2) the impact of task switching on SHD and contact dynamics. Both models show a strong linear relationship between the dynamics of SHD and contact dynamics, which exists for different initial conditions. The multiple-task-group model without task switching reveals the impacts of the number and spatial arrangements of task locations on information transmission. The task-switching model allows task-switching with a probability through contact between individuals. The model indicates that the task-switching mechanism enables a dynamical state of task-related spatial fidelity at the individual level. This spatial fidelity can assist the colony in redistributing their workforce, with consequent effects on the dynamics of spatial heterogeneity degree. The spatial fidelity of a task group is the proportion of workers who perform that task and have preferential walking styles toward their task location. Our analysis shows that the task switching rate between two tasks is an exponentially decreasing function of the spatial fidelity and contact rate. Higher spatial fidelity leads to more agents aggregating to task location, reducing contact between groups, thus making task switching more difficult. Our results provide important insights into the mechanisms that generate spatial heterogeneity and deepen our understanding of how spatial heterogeneity impacts task allocation, social interaction, and information spread.

Mathematical Modelling of Parasite Dynamics: A Stochastic Simulation-Based Approach and Parameter Estimation via Modified Sequential-Type Approximate Bayesian Computation.

The development of mathematical models for studying newly emerging and re-emerging infectious diseases has gained momentum due to global events. The gyrodactylid-fish system, like many host-parasite systems, serves as a valuable resource for ecological, evolutionary, and epidemiological investigations owing to its ease of experimental manipulation and long-term monitoring. Although this system has an existing individual-based model, it falls short in capturing information about species-specific microhabitat preferences and other biological details for different Gyrodactylus strains across diverse fish populations. This current study introduces a new individual-based stochastic simulation model that uses a hybrid τ -leaping algorithm to incorporate this essential data, enhancing our understanding of the complexity of the gyrodactylid-fish system. We compare the infection dynamics of three gyrodactylid strains across three host populations. A modified sequential-type approximate Bayesian computation (ABC) method, based on sequential Monte Carlo and sequential importance sampling, is developed. Additionally, we establish two penalised local-linear regression methods (based on L1 and L2 regularisations) for ABC post-processing analysis to fit our model using existing empirical data. With the support of experimental data and the fitted mathematical model, we address open biological questions for the first time and propose directions for future studies on the gyrodactylid-fish system. The adaptability of the mathematical model extends beyond the gyrodactylid-fish system to other host-parasite systems. Furthermore, the modified ABC methodologies provide efficient calibration for other multi-parameter models characterised by a large set of correlated or independent summary statistics.

Understanding and Quantifying Network Robustness to Stochastic Inputs.

A variety of biomedical systems are modeled by networks of deterministic differential equations with stochastic inputs. In some cases, the network output is remarkably constant despite a randomly fluctuating input. In the context of biochemistry and cell biology, chemical reaction networks and multistage processes with this property are called robust. Similarly, the notion of a forgiving drug in pharmacology is a medication that maintains therapeutic effect despite lapses in patient adherence to the prescribed regimen. What makes a network robust to stochastic noise? This question is challenging due to the many network parameters (size, topology, rate constants) and many types of noisy inputs. In this paper, we propose a summary statistic to describe the robustness of a network of linear differential equations (i.e. a first-order mass-action system). This statistic is the variance of a certain random walk passage time on the network. This statistic can be quickly computed on a modern computer, even for complex networks with thousands of nodes. Furthermore, we use this statistic to prove theorems about how certain network motifs increase robustness. Importantly, our analysis provides intuition for why a network is or is not robust to noise. We illustrate our results on thousands of randomly generated networks with a variety of stochastic inputs.

Advice to a Young Mathematical Biologist.

This paper offers advice to early-mid career researchers in Mathematical Biology from ten past and current Presidents of the Society for Mathematical Biology. The topics covered include deciding if a career in academia is right for you; finding and working with a mentor; building collaborations and working with those from other disciplines; formulating a research question; writing a paper; reviewing papers; networking; writing fellowship or grant proposals; applying for faculty positions; and preparing and giving lectures. While written with mathematical biologists in mind, it is hoped that this paper will be of use to early and mid career researchers across the mathematical, physical and life sciences, as they embark on careers in these disciplines.

Mathematical Models of Early Hepatitis B Virus Dynamics in Humanized Mice.

Analyzing the impact of the adaptive immune response during acute hepatitis B virus (HBV) infection is essential for understanding disease progression and control. Here we developed mathematical models of HBV infection which either lack terms for adaptive immune responses, or assume adaptive immune responses in the form of cytolytic immune killing, non-cytolytic immune cure, or non-cytolytic-mediated block of viral production. We validated the model that does not include immune responses against temporal serum hepatitis B DNA (sHBV) and temporal serum hepatitis B surface-antigen (HBsAg) experimental data from mice engrafted with human hepatocytes (HEP). Moreover, we validated the immune models against sHBV and HBsAg experimental data from mice engrafted with HEP and human immune system (HEP/HIS). As expected, the model that does not include adaptive immune responses matches the observed high sHBV and HBsAg concentrations in all HEP mice. By contrast, while all immune response models predict reduction in sHBV and HBsAg concentrations in HEP/HIS mice, the Akaike Information Criterion cannot discriminate between non-cytolytic cure (resulting in a class of cells refractory to reinfection) and antiviral block functions (of up to 99 % viral production 1-3 weeks following peak viral load). We can, however, reject cytolytic killing, as it can only match the sHBV and HBsAg data when we predict unrealistic levels of hepatocyte loss.

Artificial Neural Network Prediction of COVID-19 Daily Infection Count.

This study addresses COVID-19 testing as a nonlinear sampling problem, aiming to uncover the dependence of the true infection count in the population on COVID-19 testing metrics such as testing volume and positivity rates. Employing an artificial neural network, we explore the relationship among daily confirmed case counts, testing data, population statistics, and the actual daily case count. The trained artificial neural network undergoes testing in in-sample, out-of-sample, and several hypothetical scenarios. A substantial focus of this paper lies in the estimation of the daily true case count, which serves as the output set of our training process. To achieve this, we implement a regularized backcasting technique that utilize death counts and the infection fatality ratio (IFR), as the death statistics and serological surveys (providing the IFR) as more reliable COVID-19 data sources. Addressing the impact of factors such as age distribution, vaccination, and emerging variants on the IFR time series is a pivotal aspect of our analysis. We expect our study to enhance our understanding of the genuine implications of the COVID-19 pandemic, subsequently benefiting mitigation strategies.

R 0 May Not Tell Us Everything: Transient Disease Dynamics of Some SIR Models Over Patchy Environments.

This paper examines the short-term or transient dynamics of SIR infectious disease models in patch environments. We employ reactivity of an equilibrium and amplification rates, concepts from ecology, to analyze how dispersals/travels between patches, spatial heterogeneity, and other disease-related parameters impact short-term dynamics. Our findings reveal that in certain scenarios, due to the impact of spatial heterogeneity and the dispersals, the short-term disease dynamics over a patch environment may disagree with the long-term disease dynamics that is typically reflected by the basic reproduction number. Such an inconsistence can mislead the public, public healthy agencies and governments when making public health policy and decisions, and hence, these findings are of practical importance.

Modeling of Mouse Experiments Suggests that Optimal Anti-Hormonal Treatment for Breast Cancer is Diet-Dependent.

Estrogen receptor positive breast cancer is frequently treated with anti-hormonal treatment such as aromatase inhibitors (AI). Interestingly, a high body mass index has been shown to have a negative impact on AI efficacy, most likely due to disturbances in steroid metabolism and adipokine production. Here, we propose a mathematical model based on a system of ordinary differential equations to investigate the effect of high-fat diet on tumor growth. We inform the model with data from mouse experiments, where the animals are fed with high-fat or control (normal) diet. By incorporating AI treatment with drug resistance into the model and by solving optimal control problems we found differential responses for control and high-fat diet. To the best of our knowledge, this is the first attempt to model optimal anti-hormonal treatment for breast cancer in the presence of drug resistance. Our results underline the importance of considering high-fat diet and obesity as factors influencing clinical outcomes during anti-hormonal therapies in breast cancer patients.

NIAID/SMB Workshop on Multiscale Modeling of Infectious and Immune-Mediated Diseases.

On July 19th, 2023, the National Institute of Allergy and Infectious Diseases co-organized a workshop with the Society of Mathematical Biology, with the authors of this paper as the organizing committee. The workshop, "Bridging multiscale modeling and practical clinical applications in infectious diseases" sought to create an environment for mathematical modelers, statisticians, and infectious disease researchers and clinicians to exchange ideas and perspectives.

Parameter Identifiability in PDE Models of Fluorescence Recovery After Photobleaching.

Identifying unique parameters for mathematical models describing biological data can be challenging and often impossible. Parameter identifiability for partial differential equations models in cell biology is especially difficult given that many established in vivo measurements of protein dynamics average out the spatial dimensions. Here, we are motivated by recent experiments on the binding dynamics of the RNA-binding protein PTBP3 in RNP granules of frog oocytes based on fluorescence recovery after photobleaching (FRAP) measurements. FRAP is a widely-used experimental technique for probing protein dynamics in living cells, and is often modeled using simple reaction-diffusion models of the protein dynamics. We show that current methods of structural and practical parameter identifiability provide limited insights into identifiability of kinetic parameters for these PDE models and spatially-averaged FRAP data. We thus propose a pipeline for assessing parameter identifiability and for learning parameter combinations based on re-parametrization and profile likelihoods analysis. We show that this method is able to recover parameter combinations for synthetic FRAP datasets and investigate its application to real experimental data.

A Mathematical Model of Stroma-Supported Allometric Tumor Growth.

Mounting empirical research suggests that the stroma, or interface between healthy and cancerous tissue, is a critical determinate of cancer invasion. At the same time, a cancer cell's location and potential to proliferate can influence its sensitivity to cancer treatments. In this paper, we use ordinary differential equations to develop spatially structured models for solid tumors wherein the growth of tumor components is coordinated. The model tumors feature two components, a proliferating peripheral growth region, which potentially includes a mix of cancerous and noncancerous stroma cells, and a solid tumor core. Mathematical and numerical analysis are used to investigate how coordinated expansion of the tumor growth region and core can influence overall growth dynamics in a variety of tumor types. Model assumptions, which are motivated by empirical and in silico solid tumor research, are evaluated through comparison to tumor volume data and existing models of tumor growth.

Combinatorial Cooperativity in miR200-Zeb Feedback Network can Control Epithelial-Mesenchymal Transition.

Carcinomas often utilize epithelial-mesenchymal transition (EMT) programs for cancer progression and metastasis. Numerous studies report SNAIL-induced miR200/Zeb feedback circuit as crucial in regulating EMT by placing cancer cells in at least three phenotypic states, viz. epithelial (E), hybrid (h-E/M), mesenchymal (M), along the E-M phenotypic spectrum. However, a coherent molecular-level understanding of how such a tiny circuit controls carcinoma cell entrance into and residence in various states is lacking. Here, we use molecular binding data and mathematical modeling to report that the miR200/Zeb circuit can essentially utilize combinatorial cooperativity to control E-M phenotypic plasticity. We identify minimal combinatorial cooperativities that give rise to E, h-E/M, and M phenotypes. We show that disrupting a specific number of miR200 binding sites on Zeb as well as Zeb binding sites on miR200 can have phenotypic consequences-the circuit can dynamically switch between two (E, M) and three (E, h-E/M, M) phenotypes. Further, we report that in both SNAIL-induced and SNAIL knock-out miR200/Zeb circuits, cooperative transcriptional feedback on Zeb as well as Zeb translation inhibition due to miR200 are essential for the occurrence of intermediate h-E/M phenotype. Finally, we demonstrate that SNAIL can be dispensable for EMT, and in the absence of SNAIL, the transcriptional feedback can control cell state transition from E to h-E/M, to M state. Our results thus highlight molecular-level regulation of EMT in miR200/Zeb circuit and we expect these findings to be crucial to future efforts aiming to prevent EMT-facilitated dissemination of carcinomas.

Adaptation to numerosity affects the pupillary light response.

We recently showed that the gain of the pupillary light response depends on numerosity, with weaker responses to fewer items. Here we show that this effect holds when the stimuli are physically identical but are perceived as less numerous due to numerosity adaptation. Twenty-eight participants adapted to low (10 dots) or high (160 dots) numerosities and subsequently watched arrays of 10-40 dots, with variable or homogeneous dot size. Luminance was constant across all stimuli. Pupil size was measured with passive viewing, and the effects of adaptation were checked in a separate psychophysical session. We found that perceived numerosity was systematically lower, and pupillary light responses correspondingly smaller, following adaptation to high rather than low numerosities. This is consistent with numerosity being a primary visual feature, spontaneously encoded even when task irrelevant, and affecting automatic and unconscious behaviours like the pupillary light response.

Enumeration of Rooted Binary Unlabeled Galled Trees.

Rooted binary galled trees generalize rooted binary trees to allow a restricted class of cycles, known as galls. We build upon the Wedderburn-Etherington enumeration of rooted binary unlabeled trees with n leaves to enumerate rooted binary unlabeled galled trees with n leaves, also enumerating rooted binary unlabeled galled trees with n leaves and g galls, 0 ⩽ g ⩽ ⌊ n - 1 2 ⌋ . The enumerations rely on a recursive decomposition that considers subtrees descended from the nodes of a gall, adopting a restriction on galls that amounts to considering only the rooted binary normal unlabeled galled trees in our enumeration. We write an implicit expression for the generating function encoding the numbers of trees for all n. We show that the number of rooted binary unlabeled galled trees grows with 0.0779 ( 4 . 8230 n ) n - 3 2 , exceeding the growth 0.3188 ( 2 . 4833 n ) n - 3 2 of the number of rooted binary unlabeled trees without galls. However, the growth of the number of galled trees with only one gall has the same exponential order 2.4833 as the number with no galls, exceeding it only in the subexponential term, 0.3910 n 1 2 compared to 0.3188 n - 3 2 . For a fixed number of leaves n, the number of galls g that produces the largest number of rooted binary unlabeled galled trees lies intermediate between the minimum of g = 0 and the maximum of g = ⌊ n - 1 2 ⌋ . We discuss implications in mathematical phylogenetics.

Discretised Flux Balance Analysis for Reaction-Diffusion Simulation of Single-Cell Metabolism.

Metabolites have to diffuse within the sub-cellular compartments they occupy to specific locations where enzymes are, so reactions could occur. Conventional flux balance analysis (FBA), a method based on linear programming that is commonly used to model metabolism, implicitly assumes that all enzymatic reactions are not diffusion-limited though that may not always be the case. In this work, we have developed a spatial method that implements FBA on a grid-based system, to enable the exploration of diffusion effects on metabolism. Specifically, the method discretises a living cell into a two-dimensional grid, represents the metabolic reactions in each grid element as well as the diffusion of metabolites to and from neighbouring elements, and simulates the system as a single linear programming problem. We varied the number of rows and columns in the grid to simulate different cell shapes, and the method was able to capture diffusion effects at different shapes. We then used the method to simulate heterogeneous enzyme distribution, which suggested a theoretical effect on variability at the population level. We propose the use of this method, and its future extensions, to explore how spatiotemporal organisation of sub-cellular compartments and the molecules within could affect cell behaviour.

Second-Order Effects of Chemotherapy Pharmacodynamics and Pharmacokinetics on Tumor Regression and Cachexia.

Drug dose response curves are ubiquitous in cancer biology, but these curves are often used to measure differential response in first-order effects: the effectiveness of increasing the cumulative dose delivered. In contrast, second-order effects (the variance of drug dose) are often ignored. Knowledge of second-order effects may improve the design of chemotherapy scheduling protocols, leading to improvements in tumor response without changing the total dose delivered. By considering treatment schedules with identical cumulative dose delivered, we characterize differential treatment outcomes resulting from high variance schedules (e.g. high dose, low dose) and low variance schedules (constant dose). We extend a previous framework used to quantify second-order effects, known as antifragility theory, to investigate the role of drug pharmacokinetics. Using a simple one-compartment model, we find that high variance schedules are effective for a wide range of cumulative dose values. Next, using a mouse-parameterized two-compartment model of 5-fluorouracil, we show that schedule viability depends on initial tumor volume. Finally, we illustrate the trade-off between tumor response and lean mass preservation. Mathematical modeling indicates that high variance dose schedules provide a potential path forward in mitigating the risk of chemotherapy-associated cachexia by preserving lean mass without sacrificing tumor response.