Mathematics

In many natural dynamical systems, tipping points emerge where slowly changing external factors induce shifts from one state to another, known as critical transitions or bifurcations. In epidemic systems, detecting these transitions is important for outbreak preparedness. Various statistical indicators like lag-1 autocorrelation and variance have been employed as early warning signals (EWS) of critical transitions by showing how the dynamics of the system exhibit critical slowing down (CSD) near the transition. More recently, deep learning-based algorithms have been developed that attempt to more accurately detect and classify bifurcations; however, they are limited by the systems simulated in the training dataset, the length of the time series, and have not been rigorously evaluated on real-world epidemic data spanning multiple geographic regions. In this paper, we evaluate existing deep learning-based bifurcation prediction models, on a CDC-curated dataset containing weekly influenza hospital admissions across all fifty U.S. states and the nation as a whole. We develop a pipeline for data preparation, reproduction number (Rt) estimation using the EpyEstim package (a Python implementation of the EpiEstim framework [4]), and bifurcation point identification. We systematically evaluate the models under several windowing strategies, including rolling, expanding, and fixed-length windows, and report accuracy, precision, recall, F1 score, and specificity. Our best-performing configuration (expanding window, length 100) achieves an accuracy of 90.09%, precision of 72.86%, recall of 96.23%, and F1 score of 82.93%. We additionally discuss Markov Regime Switching models as an alternative framework for detecting epidemic transitions. Our findings indicate that these models can generalize to real-world, state-level influenza data and are capable of online/real-time prediction.

In this study, we develop a mathematical model of the hypothalamic pituitary adrenal (HPA) axis and systematically investigate the model complexity required to reproduce observed hormone dynamics. The model is fitted to ACTH and cortisol measurements from five healthy volunteers reported in the literature. We show that simultaneous reproduction of both hormone profiles requires a delay differential equation model with at least five state variables among the model architectures considered. Following fitting the volunteer data, we investigate the response of the model to sustained stress using bifurcation analysis. The coupled system exhibits bifurcations over a range of stress intensities. This analysis reveals bifurcations whose critical stress thresholds depend strongly on the modulatory state of the system. These results provide a mathematical framework for investigating transitions between physiological and dysregulated HPA axis dynamics, which may be related to states of health and disease.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Mitochondrial reactive oxygen species (ROS) are implicated in cardiac dysfunction, but complex dynamic interactions between ROS, intracellular calcium, and electrophysiology make it difficult to resolve mechanisms experimentally. We developed the Zukowski excitation-contraction-mitochondrial-ROS (ECM-ROS) model, a computational model that couples human ventricular electrophysiology with mitochondrial calcium handling, energetics, ROS production, and ROS dependent modulation of the ryanodine receptor, SERCA, L-type calcium current, and late sodium current. Model predictions were evaluated against independent experimental measurements of ROS induced changes in action potential duration and intracellular calcium. We then used the model to characterize ROS calcium feedback across physiological pacing rates and applied it to doxorubicin induced oxidative stress. Increasing cytosolic calcium produced nonlinear amplification of mitochondrial activity and cytosolic ROS, with greater amplification at faster pacing rates. Moderate ROS elevation initially enhanced calcium release, whereas greater oxidative stress depleted sarcoplasmic reticulum calcium stores. During simulated doxorubicin exposure, faster pacing lowered the level of oxidative stress required to produce calcium dysregulation and action potential prolongation. Mitochondrial ROS and calcium accumulation preceded detectable electrophysiological remodeling. In populations incorporating variability in ion channel expression, neither doxorubicin induced oxidative stress nor subtle pharmacological IKr block alone produced early afterdepolarizations. However, early afterdepolarizations were predicted to emerge with combined perturbations. These findings establish a computational framework for investigating bidirectional coupling among mitochondrial function, ROS, calcium handling, and human cardiac electrophysiology. The model generates testable predictions regarding the early detection of oxidative cardiac injury and the physiological conditions that increase susceptibility to arrhythmia.

To reduce health risks for human space exploration, it is important to model the effects of deep-space radiation on biological systems. The budding yeast Saccharomyces cerevisiae is a common model organism in space radiobiology, but little is known about how radiation damage to individual cells translates into the population-level effects that experiments measure. The Agent-Based Model for Microbial Populations Exposed to Radiation (AMMPER) is a spatially explicit agent-based model that simulates the effects of ionizing radiation on S. cerevisiae cells and produces population-level results. Its goal is to improve experiment design and data interpretation. AMMPER-1 demonstrated that agent-based simulations can qualitatively reproduce the effects of proton radiation on the growth rate of wild-type and DNA-repair mutant yeast. Here, we present AMMPER-2, which addresses limitations of the first model and expands its capabilities. AMMPER-2 incorporates the dynamics of the redox dye alamarBlue, a data type commonly generated in spaceflight microbiology experiments, and models the fate of reactive oxygen species in more detail. It also adds an improved graphical interface and runs faster. AMMPER-2 remains easy to execute in Python on a standard personal laptop, providing space biologists a tool to explore experimental conditions, generate hypotheses, and interpret data from microbial space radiation experiments.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [≥] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

We propose a stochastic epidemic model to describe the spread of avian influenza in a network of poultry farms. We consider three main paths of infection transmission between poultry farms: from nearby farms, mainly through airborne diffusion; from farms belonging to the same company, through shared veterinarians, forage providers and similar; from undetected small farms or wildlife. From the parameter estimation, based on a modified Expectation-Maximization algorithm, we infer that approximately 63% of the farms were infected from nearby ones, with an infection force declining with distance; of these, more than one third belonged to the same company of the estimated infector. About 20% were infected from premises belonging to the same company but farther away than the distance threshold of 2 km; the remaining ones from wildlife or unidentified sources. These estimates have been validated by a comparison with genetic data, available for a subset of the farms: the genetic distance between two farms identified, with high probability, as an infector-infectee pair is much lower than between random pairs. These results may help in the implementation of tailored prevention and control measures for future outbreaks.

In this work, we investigate the homogeneity of performance in binary classification between the Rotterdam phenotypes and phenotype prediction. For this purpose, a multi-task neural network architecture is designed to improve the prediction accuracy of phenotypes while maintaining the accuracy of binary classification. Using a clinical PCOS dataset which contains unlabeled samples of patients that cannot be attributed to any of the known phenotypes, three models are evaluated: optimized binary Random Forest (RF), six-class weighted class RF, and multi-task neural network (with a common hidden layer and two heads for binary and six-class classifications) optimized by a genetic algorithm. Techniques used for evaluation include ten-fold out-of-fold predictions, paired McNemar tests, bootstrap tests, ablation study, and clustering (Adjusted Rand Index, ARI). Binary RF model reached 90.2% accuracy; however, recall of different phenotypes was significantly unbalanced - from 75.0% of phenotype B to 97.8% of phenotype A. This difference of 22.8 % points was observed (23.3 points) even after excluding phenotype-specific features. While direct phenotype prediction improved recall for common phenotypes, it decreased recall for phenotype B to 43.8% and unclassified cases to 31.5%. Our multi-task model preserved 90.2% binary accuracy and significantly recovered phenotype B recall to 87.5% and unclassified cases to 79.6% (p<0.05). The cluster analysis confirmed the weak correlation between the results and labels (ARI=0.045). The ablation analysis demonstrates that feature selection and class weighting drive this improvement, isolating the specific components responsible for performance gains. We establish that aggregate binary accuracy does not reflect the heterogeneous detection of PCOS phenotypes. Direct phenotype prediction alone redistributes the classification errors, but the multi-task architecture successfully improves the detection of rare phenotypes without compromising binary accuracy. Phenotype-disaggregated evaluations, supported by significance testing and ablation, are essential for clinical classifier research.

Simplified spiking neuron models such as the Adaptive Exponential Integrate-and-Fire (AdEx) model reproduce the firing behaviour of diverse neuron types at a fraction of the cost of biophysically detailed models, making them attractive for large-scale brain simulation. Fitting their parameters to recordings, however, still relies on derivative-free search such as grid search or evolutionary algorithms, because the discrete spike-and-reset mechanism renders these models non-differentiable. Surrogate gradients, used in deep spiking networks, replace the spike derivative with a smooth approximation on the backward pass. Whether they also enable efficient gradient-based parameter estimation for single simplified neurons is untested. We implement a differentiable AdEx model with a surrogate-gradient variant in the Jaxley framework. The implementation is on par with target-compiled C++ Brian2 code while remaining end-to-end differentiable, vmap-batchable and GPU/TPU-portable. It is 200 times faster than Brian2's default runtime mode. Across a recovery-radius benchmark of 2,700 paired runs, gradient-based fitting never matched a gradient-free Nelder-Mead baseline. At sub- threshold levels, gradient descent instead beats Nelder-Mead. The obstacle is therefore the spike mechanism and the loss geometry it induces, not gradient descent itself.

Traditional cell tracking methods focus on fully automated tracking. The quality of the generated tracks, however, is heavily influenced by factors such as cell speed and density. In complex scenarios where track correctness is essential, an efficient and reliable tool for manual track correction remains a key component of the cell tracking process. Most tracking tools that allow for such manual correction do so only as a post-processing step, using human feedback to correct one track error at a time -- without improving the quality of other tracks in the dataset. This highlights the need for a more efficient track correction method that iteratively integrates human feedback into the automated tracking process. In this paper we propose a human-in-the-loop approach to cell tracking and investigate the extent to which a tracking algorithm can learn from human input during the tracking process. We combine a graph representation of cell detections, Dijkstra's shortest path algorithm, and a neural network that predicts cell displacement to generate initial track suggestions. We introduce a tool that enables human feedback on individual tracks suggested by the algorithm; this feedback is then propagated to improve subsequent track generation. Compared to the traditional approach of fully automated tracking with correction as a post-processing step, our human-in-the-loop tracking framework reduces the number of corrections required. Finally, we show how this approach can be integrated with existing state-of-the-art tracking algorithms by using pre-generated tracks as a prior, improving the quality of tracks suggested to the human annotator.

In recent years, the issue of reproducibility in scientific studies has been at the forefront of scientific discussion. Many of these issues arise due to poor employment of statistical methods, making it crucial to address these knowledge gaps for frequently used methods. Elastic net regression (ENET) is one such algorithm used frequently in biology and public health due to its ability to select important predictors (also known as "features") from large data sets. However, ENET struggles with consistent feature selection and hyperparameter choice in these applications due to low observation counts and noisy data. This limitation hinders ENET's effectiveness in biological research due to low model reproducibility and predictor interpretability. We propose the Instability Quotient Criterion (IQC) as a metric to quantify the consistency of selected features for feature selection models such as ENET, with a specific emphasis on biological applications. IQC leverages ensemble modeling and Shannon's entropy in order to determine the stability of selected features derived from ENET models. We demonstrate IQC's robustness via simulation studies, showing that IQC is capable of detecting significant decreases in model stability with lower observation counts, higher dimensionality, and higher proportions of noisy signal. We additionally utilize simulation to provide a naive threshold for low versus high stability. Finally, we applied IQC to a case study to illustrate how our method can be integrated into a larger analysis pipeline when looking at microarray gene expression data. IQC indicated stability issues for differentially-expressed genes when classifying acute leukemia subtypes in spite of high classification success. Lastly, we discuss how model instability impedes interpretation of biologically significant genes for discerning cancer subtypes. Overall, we present IQC as a robust evaluation metric that can easily be incorporated into high-dimensional biological analyses.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Quantitative characterization of cell proliferation is central to preclinical drug discovery. Here, we evaluated Random Forest (RF) regression for predicting confluence-based cell growth trends using data from human cancer cell lines and benchmarked its performance against the widely used logistic and Gompertz models. The RF model achieved higher predictive accuracy within the evaluated dataset, particularly at higher initial seeding densities. We further extended the RF framework to predict growth curves across seeding densities. Collectively, these findings highlight the potential of data-driven modeling to complement conventional mathematical approaches and support more informed cell-culture design.

Dynamic systems often evolve faster than reliable measurements can be collected. More frequent indirect observations may fill these gaps, but their noise and differing scales complicate estimation of the underlying state. We developed Bayesian Anchored Latent Learning (BALL) to combine sparse reference measurements with denser records on an interpretable scale. A bidirectional transformer teacher reconstructs trajectories from complete training histories, and a forward-only transformer student learns from those reconstructions and observed measurements. Across 30 simulations, BALL reduced composite trajectory error by 7% relative to the same transformer trained directly from measurements, achieving lower error across all datasets. In an applied example using psychiatric records from 1,609 patients, BALL improved the prediction of withheld questionnaire scores compared with direct transformer training at 21- and 28-day measurement schedules. Requests guided by estimated uncertainty further improved prediction. BALL connects retrospective learning with prospective state estimation and provides a framework for directing additional measurement.

Allergic conjunctivitis is an inflammation of the conjunctiva caused by allergen; it is common disorder affecting up to 40% of the population. In this work, we study seasonal allergic conjunctivitis (SAC), also called "hay fever eyes", which is caused by exposure to airborne pollens. We develop a mathematical model quantifying the ocular immune system response to the allergens. First, we present a simplified qualitative description of the immunopathogenesis of SAC. Then, we express each chosen immunopathological mechanism mathematically to construct a system of thirty-one ordinary differential equations. We compare summary statistics of the predicted observable immune signals to experimental measurements and find our model captures key qualitative features of SAC progression. We then compare our predicted time series of histamine concentration to symptom scores and find a strong correlation suggesting the model predicts relevant clinically trends. Next, we calibrate the model through multi-step process. We find the most influential parameters are the production and depletion rates of IL-4, and the production rates of IL-5 and IL-8. These cytokines are targeted in treatments for asthma, atopic dermatitis, and severe eosinophilic associated disorder and suggest potential therapeutic targets for SAC. Our calibrated model mimics most of the summary statistics of the experimentally observable immune signals with discrepancies for IL-5 and IL-13 indicating that additional immunopathological mechanisms could be important.

Computer vision approaches to disease motor assessment are limited by clinical video dataset scarcity and distributional imbalance, motivating interest in generative video models as a potential source of training data. We introduce a three-component framework for evaluating the clinical fidelity of synthetic medical motion data, assessing visual fidelity, biomechanical stability, and pathophysiological accuracy, and apply it to videos generated by the text-to-video model Sora 2 depicting Parkinsonian hand-fisting motions across MDSUPDRS severity levels 0-2 (normal to mild-moderate motor impairment). We further use real-versus-real bootstrapping to compare synthetic-real differences against the natural variability of real clinical videos. Although generated videos appear visually realistic, quantitative analysis shows substantial divergence from real clinical recordings, including weaker severity-dependent patterns of bradykinesia and hypokinesia, less realistic tremor rhythms, and substantial temporal instability. These findings reveal a motion-domain gap between Sora-2-generated and real clinical motor behavior in this Parkinsonian hand-motion setting and demonstrate that visual realism alone is insufficient to establish the clinical utility of synthetic video data for training or evaluating medical computer vision systems.

Nitrogen dioxide (NO2) is a widespread air pollutant associated with the risk and prevalence of various respiratory infectious diseases. Epidemiological studies have reported associations between NO2 exposure and COVID-19 incidence, and biological evidence indicates that NO2 exposure increases vulnerability to infection. Motivated by these findings and by the scarcity of mechanistic epidemic models that incorporate air pollutant data, here we introduce an NO2- and mobility-informed metapopulation model and study the first COVID-19 wave in Oklahoma. County-level NO2 concentrations enter the infection rates through a first-order parameterization, alongside high-resolution mobility patterns capturing heterogeneous inter-county connectivity. Using numerical simulations, we demonstrate that higher environmental sensitivity produces earlier and more synchronized outbreaks across counties, while lower sensitivity yields delayed, asynchronous trajectories. When fitted to county-level COVID-19 incidence data, our model outperformed the classic SIR model (fitted independently in each county) in over 75% of counties in Oklahoma across all standard error metrics. The improvements were even more pronounced in forecasting tasks across epidemic stages, with county-level win rates against the SIR model of 59.74% (surge), 87.01% (peak), 80.52% (post-peak), and 88.31% (decline) by nRMSE. To isolate the contribution of NO2, we compared our model against a mobility-only model with NO2 removed and homogeneous infection rates across counties. The NO2 term improved county-level forecasts during the surge (94.81% of counties by nRMSE), post-peak (88.31%), and peak (62.34%) stages, though not during the early and decline stages, suggesting that NO2 information can generate better forecasts than mobility alone during periods of high incidence. Our results demonstrate that environmental exposure can enter the transmission structure of a mechanistic epidemic model directly and that doing so yields measurable forecasting gains.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Simulation models are widely used to support counterfactual inference in settings where direct experimental validation is difficult or infeasible. However, for phenomenological models calibrated to observational time series, goodness-of-fit does not necessarily establish that the structural assumptions underlying such inferences are uniquely supported by the data. Here, we examine this issue using a published COVID-19 transmission model for Japan that incorporates empirical inputs including human mobility, variant prevalence, and vaccination and has been used to estimate the effects of alternative vaccination scenarios. Using the publicly archived code together with a reconstructed substitute for a restricted-access input file required by the released implementation, we first reproduced the published model behaviour. We then performed two complementary computational stress tests. First, the vaccination-related input was suppressed and the remaining parameters were re-estimated against the original epidemic trajectory. Second, as a complementary test of trajectory-level flexibility, the model was refitted to synthetic target series with qualitatively different temporal structures while retaining the empirical mobility and variant inputs. The model retained substantial fitting capacity in both tests: suppression of the vaccination-related input caused only a limited deterioration in fit to the original trajectory, while qualitatively different synthetic trajectories could also be reproduced through parameter adjustment. Variation in fitted parameters across these conditions was also consistent with compensatory parameter behaviour. These results do not imply that vaccination had no effect or that the original model is invalid. Rather, they provide empirical evidence that, in this case, goodness-of-fit to the observed epidemic trajectory alone does not uniquely support the structural interpretation required for vaccination-related counterfactual inference. The study illustrates the value of reproducible computational stress testing as a complement to, rather than a substitute for, formal identifiability or model-comparison analyses.