Predator-prey interactions are one of ecology's central research themes, but with many interdisciplinary implications across the social and natural sciences. Here we consider an often-overlooked species in these interactions, namely parasites. We first show that a simple predator-prey-parasite model, inspired by the classical Lotka-Volterra equations, fails to produce a stable coexistence of all three species, thus failing to provide a biologically realistic outcome. To improve this, we introduce free space as a relevant eco-evolutionary component in a new mathematical model that uses a game-theoretical payoff matrix to describe a more realistic setup. We then show that the consideration of free space stabilizes the dynamics by means of cyclic dominance that emerges between the three species. We determine the parameter regions of coexistence as well as the types of bifurcations leading to it by means of analytical derivations as well as by means of numerical simulations. We conclude that the consideration of free space as a finite resource reveals the limits of biodiversity in predator-prey-parasite interactions, and it may also help us in the determination of factors that promote a healthy biota.

This work presents a mathematical model, based on partial differential equations, that analyzes the inflammatory stage of atherosclerosis. Four leading players are taken into consideration: Low Density Lipoproteins, oxidized Low Density Lipoproteins, immune cells and the inflammatory cytokines. In addition to this, the permeability of the endothelial layer is taken into account in the model. A stability analysis of the fixed points of the kinetic system is presented in details followed by the proof of existence of traveling wave solutions of the system of partial differential equations. The mathematical analysis leads to a biological interpretation. We distinguish three main cases of the disease state that correlate with the permeability of the endothelial layer. In fact, having a low permeability indicates the disease free state since no chronic inflammatory reaction occurs due to the non initiation of the inflammation. With intermediate permeability, a wave propagation corresponding to a chronic inflammatory reaction might occur whether the initial perturbation overcomes a threshold or not. With high permeability, even a small perturbation of the disease free state leads to a chronic inflammatory reaction represented by a wave propagation. We perform numerical simulations of the solutions to illustrate the biological results.

Mesenchymal cell motility is mainly regulated by two members of the Rho-family of GTPases, called Rac and Rho. The mutual inhibition exerted by these two proteins on each other's activation and the promotion of Rac activation by an adaptor protein called paxillin have been implicated in driving cellular polarization comprised of front (high active Rac) and back (high active Rho) during cell migration. Mathematical modeling of this regulatory network has previously shown that bistability is responsible for generating a spatiotemporal pattern underscoring cellular polarity called wave-pinning when diffusion is included. We previously developed a 6V reaction-diffusion model of this network to decipher the role of Rac, Rho and paxillin (along with other auxiliary proteins) in generating wave-pinning. In this study, we simplify this model through a series of steps into an excitable 3V ODE model comprised of one fast variable (the scaled concentration of active Rac), one slow variable (the maximum paxillin phosphorylation rate - turned into a variable) and a very slow variable (a recovery rate - also turned into a variable). We then explore, through slow-fast analysis, how excitability is manifested by showing that the model can exhibit relaxation oscillations (ROs) as well as mixed-mode oscillations (MMOs) whose underlying dynamics are consistent with a delayed Hopf bifurcation with a canard explosion. By reintroducing diffusion and the scaled concentration of inactive Rac into the model, we obtain a 4V PDE model that generates several unique spatiotemporal patterns that are relevant to cell motility. These patterns are then characterized and their impact on cell motility are explored by employing the cellular potts model (CPM). Our results reveal that wave pinning produces purely very directed motion in CPM, while MMOs allow for meandering and non-motile behaviors to occur. This highlights the role of MMOs as a potential mechanism for mesenchymal cell motility.

Chronic obstructive pulmonary disease (COPD) is a highly prevalent lung disease characterized by chronic inflammation and tissue remodeling possibly induced by unusual interactions between fibrocytes and CD8+ T lymphocytes in the peribronchial area. To investigate this phenomenon, we developed a probabilistic cellular automata type model where the two types of cells follow simple local interaction rules taking into account cell death, proliferation, migration and infiltration. We conducted a rigorous mathematical analysis using multiscale experimental data obtained in control and disease conditions to estimate the model's parameters accurately. The simulation of the model is straightforward to implement, and two distinct patterns emerged that we can analyse quantitatively. In particular, we show that the change in fibrocyte density in the COPD condition is mainly the consequence of their infiltration into the lung during exacerbations, suggesting possible explanations for experimental observations in normal and COPD tissue. Our integrated approach that combines a probabilistic cellular automata model and experimental findings will provide further insights into COPD in future studies.

Cell-based models provide a helpful approach for simulating complex systems that exhibit adaptive, resilient qualities, such as cancer. Their focus on individual cell interactions makes them a particularly appropriate strategy to study cancer therapies' effects, which are often designed to disrupt single-cell dynamics. In this work, we propose them as viable methods for studying the time evolution of cancer imaging biomarkers (IBM). We propose a cellular automata model for tumor growth and three different therapies: chemotherapy, radiotherapy, and immunotherapy, following well-established modeling procedures documented in the literature. The model generates a sequence of tumor images, from which a time series of two biomarkers: entropy and fractal dimension, is obtained. Our model shows that the fractal dimension increased faster at the onset of cancer cell dissemination. At the same time, entropy was more responsive to changes induced in the tumor by the different therapy modalities. These observations suggest that the prognostic value of the proposed biomarkers could vary considerably with time. Thus, it is essential to assess their use at different stages of cancer and for different imaging modalities. Another observation derived from the results was that both biomarkers varied slowly when the applied therapy attacked cancer cells scattered along the automatons' area, leaving multiple independent clusters of cells at the end of the treatment. Thus, patterns of change of simulated biomarkers time series could reflect on essential qualities of the spatial action of a given cancer intervention.

Within-host SARS-CoV-2 modelling studies have been published throughout the COVID-19 pandemic. These studies contain highly variable numbers of individuals and capture varying timescales of pathogen dynamics; some studies capture the time of disease onset, the peak viral load and subsequent heterogeneity in clearance dynamics across individuals, while others capture late-time post-peak dynamics. In this study, we curate multiple previously published SARS-CoV-2 viral load data sets, fit these data with a consistent modelling approach, and estimate the variability of in-host parameters including the basic reproduction number, R0, as well as the best-fit eclipse phase profile. We find that fitted dynamics can be highly variable across data sets, and highly variable within data sets, particularly when key components of the dynamic trajectories (e.g. peak viral load) are not represented in the data. Further, we investigated the role of the eclipse phase time distribution in fitting SARS-CoV-2 viral load data. By varying the shape parameter of an Erlang distribution, we demonstrate that models with either no eclipse phase, or with an exponentially-distributed eclipse phase, offer significantly worse fits to these data, whereas models with less dispersion around the mean eclipse time (shape parameter two or more) offered the best fits to the available data across all data sets used in this work. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

The intracranial pressure (ICP) curve with its different peaks has been comprehensively studied, but the exact physiological mechanisms behind its morphology has not been revealed. If the pathophysiology behind deviations from the normal ICP curve form could be identified, it could be vital information to diagnose and treat each single patient. A mathematical model of the hydrodynamics in the intracranial cavity over single heart cycles was developed. A Windkessel model approach was generalized but the unsteady Bernoulli equation was utilized for blood flow and CSF flow. This is a modification of earlier models using the extended and simplified classical Windkessel analogies to a model that is based on mechanisms rooted in the laws of physics. The improved model was calibrated with patient data for cerebral arterial inflow, venous outflow, cerebrospinal fluid (CSF), and ICP over one heart cycle from 10 neuro-intensive care unit patients. A priori model parameter values were obtained by considering patient data and values taken from earlier studies. These values were used as an initial guess for an iterated constrained-ODE (ordinary differential equation) optimization problem with cerebral arterial inflow data as input into the system of ODEs. The optimization routine found patient-specific model parameter values that produced model ICP curves that showed excellent agreement with clinical measurements while model venous and CSF flow were within a physiologically acceptable range. The improved model and the automated optimization routine gave better model calibration results compared to previous studies. Moreover, patient-specific values of physiologically important parameters like intracranial compliance, arterial and venous elastance, and venous outflow resistance were determined. The model was used to simulate intracranial hydrodynamics and to explain the underlying mechanisms of the ICP curve morphology. Sensitivity analysis showed that the order of the three main peaks of the ICP curve was affected by a decrease in arterial elastance, a large increase in resistance to arteriovenous flow, an increase in venous elastance, or a decrease in resistance to CSF flow in the foramen magnum; and the frequency of oscillations were notably affected by intracranial elastance. In particular, certain pathological peak patterns were caused by these changes in physiological parameters. To the best of our knowledge, there are no other mechanism-based models associating the pathological peak patterns to variation of the physiological parameters.

Local environmental interactions are a major factor in determining the success of a new mutant in structured populations. Spatial variations in the concentration of genotype-specific resources change the fitness of competing strategies locally and thus can drastically change the outcome of evolutionary processes in unintuitive ways. The question is how such local environmental variations in network population structures change the condition for selection and fixation probability of an advantageous (or deleterious) mutant. We consider linear graph structures and focus on the case where resources have a spatial periodic pattern. This is the simplest model with two parameters, length scale and fitness scales, representing heterogeneity. We calculate fixation probability and fixation times for a constant population birth-death process as fitness heterogeneity and period vary. Fixation probability is affected by not only the level of fitness heterogeneity but also spatial scale of resources variations set by period of distribution T. We identify conditions for which a previously a deleterious mutant (in a uniform environment) becomes beneficial as fitness heterogeneity is increased. We observe cases where the fixation probability of both mutant and resident types are more than their neutral value, 1/N, simultaneously. This coincides with exponential increase in time to fixation which points to potential coexistence of resident and mutant types. Finally, we discuss the effect of the 'fitness shift' where the fitness function of two types has a phase difference. We observe significant increases (or decreases) in the fixation probability of the mutant as a result of such phase shift.

Understanding the mechanistic dynamics of transmission is key to designing more targeted and effective interventions to limit the spread of infectious diseases. A well-described within-host model allows explicit simulation of how infectiousness changes over time at an individual level. This can then be coupled with dose-response models to investigate the impact of timing on transmission. We collected and compared a range of within-host models used in previous studies and identified a minimally-complex model that provides suitable within-host dynamics while keeping a reduced number of parameters to allow inference and limit unidentifiability issues. Furthermore, non-dimensionalised models were developed to further overcome the uncertainty in estimates of the size of the susceptible cell population, a common problem in many of these approaches. We will discuss these models, and their fit to data from the human challenge study (see Killingley et al. (2022)) for SARS-CoV-2 and the model selection results, which has been performed using ABC-SMC. The parameter posteriors have then used to simulate viral-load based infectiousness profiles via a range of dose-response models, which illustrate the large variability of the periods of infection window observed for COVID-19.

Humans and other primates exhibit pro-social preferences for fairness. These preferences are thought to be reinforced by strong reciprocity, a policy that rewards fair actors and punishes unfair ones. Theories of fairness based on strong reciprocity have been criticized for overlooking the importance of individual differences in socially heterogeneous populations. Here, we explore the evolution of fairness in a heterogeneous population. We analyse the Ultimatum Game in cases where players' roles in the game are determined by their status. Importantly, our model allows for non-random pairing of players, and so we also explore the role played by kin selection in shaping fairness. Our kin-selection model shows that, when individuals condition their behaviour on their role in the game, fairness can be understood as either altruistic or spiteful. Altruistic fairness directs resources from less valuable members of a genetic lineage to more valuable members of the same lineage, whereas spiteful fairness keeps resources away from the competitors of the actor's high-value relatives. When individuals express fairness unconditionally it can be understood as altruistic or selfish. When it is altruistic, unconditional fairness again serves to direct resources to high-value members of genetic lineages. When it is selfish, unconditional fairness simply improves an individual's own standing. Overall, we expand kin-selection based explanations for fairness to include motivations other than spite. We show, therefore, that one need not invoke strong reciprocity to explain the advantage of fairness in heterogeneous populations.

Understanding the role of natural selection in driving evolutionary change requires accurate estimates of the strength of selection acting at the genetic level in the wild. This is challenging to achieve but may be easier in the case of populations in migration-selection balance. When two populations are at equilibrium under migration-selection balance, there exist loci whose alleles are selected different ways in the two populations. Such loci can be identified from genome sequencing by their high values of FST. This raises the question of what is the strength of selection on locally-adaptive alleles. To answer this question we analyse a 1-locus 2-allele model of a population distributed between two niches. We show by simulation of selected cases that the outputs from finite-population models are essentially the same as those from deterministic infinite-population models. We then derive theory for the infinite-population model showing the dependence of selection coefficients on equilibrium allele frequencies, migration rates, dominance and relative population sizes in the two niches. An Excel spreadsheet is provided for the calculation of selection coefficients and their approximate standard errors from observed values of population parameters. We illustrate our results with a worked example, with graphs showing the dependence of selection coefficients on equilibrium allele frequencies, and graphs showing how FST depends on the selection coefficients acting on the alleles at a locus. Given the extent of recent progress in ecological genomics, we hope our methods may help those studying migration-selection balance to quantify the advantages conferred by adaptive genes.

COVID-19, induced by the SARS-CoV-2 infection, has caused an unprecedented pandemic in the world. New variants of the virus have emerged and dominated the virus population. In this paper, we develop a multi-strain model with asymptomatic transmission to study how the asymptomatic or pre-symptomatic infection influences the transmission between different strains and control strategies that aim to mitigate the pandemic. Both analytical and numerical results reveal that the competitive exclusion principle still holds for the model with the asymptomatic transmission. By fitting the model to the COVID-19 case and viral variant data in the US, we show that the omicron variants are more transmissible but less fatal than the previously circulating variants. The basic reproduction number for the omicron variants is estimated to be 11.15, larger than that for the previous variants. Using mask mandate as an example of non-pharmaceutical interventions, we show that implementing it before the prevalence peak can significantly lower and postpone the peak. The time of lifting the mask mandate can affect the emergence and frequency of subsequent waves. Lifting before the peak will result in an earlier and much higher subsequent wave. Caution should also be taken to lift the restriction when a large portion of the population remains susceptible. The methods and results obtained her e may be applied to the study of the dynamics of other infectious diseases with asymptomatic transmission using other control measures.

Estimating microbial mutation rates is an essential task in evolutionary biology, with wide range applications in related fields such as virology, epidemiology, clinic and public health, and antibiotic research. Significant progress has been made on this research since 1943 when Luria-Delbrück fluctuation analysis was first introduced. However, existing estimators of mutation rates are heavily reliant on model assumptions in fluctuation analysis, and become less applicable to real microbial experiments which deviate from the model assumptions. To overcome this difficulty, we propose to model fluctuation experimental data by a two-type Markov branching process (MBP) and use approximate Bayesian computation (ABC) to estimate the mutation probability parameters. Such an ABC-based mutation rate estimator is based on intensive simulations from the mutation process, thereby taking advantage of modern computing power. Most importantly, its likelihood-free feature allows more complex and realistic setups of the mutation process, especially when the distribution of the number of mutants cannot be easily derived. To further improve computation efficiency, we use a Gaussian process surrogate to substitute the simulator in the ABC algorithm, and call the resulting estimator GPS-ABC. Simulation studies show that, when used to estimate constant mutation rate in MBP, ABC-based estimators generally outperform traditional moment or likelihood-based estimators. When mutations occur in two stages, i.e., in MBP with a piece-wise constant mutation rate function, traditional mutation rate estimators become not applicable, yet GPS-ABC still achieves reasonable estimates. Finally, the proposed GPS-ABC estimator is used to analyze real fluctuation experimental datasets for studying drug resistance.

In a long-lived organism with a modular architecture, such as trees, somatic mutations accumulate throughout the long lifespan and result in genetic mosaicism in each module within the same individual. In recent years, next-generation sequencing technology has provided a snapshot of such intra-organismal genetic variability. However, the dynamic processes underlying the accumulation and expansion of somatic mutations during the growth remain poorly understood. In this study, we constructed a model to describe these processes in a form that can be applied to a real tree. Given that the proliferation dynamics of meristematic cells vary across plant species, multiple possible processes for elongation and branching were comprehensively expressed in our model. Using published data from a poplar tree, we compared the prediction of the models with the observation and explained the cell lineage dynamics underlying somatic mutations accumulation that were not evident from the snapshot of the sequenced data. We showed that the somatic genetic drift during growth increases inter-meristem mosaicism, resulting in genetically distinct branches and less integrity within an individual tree. We also showed that the somatic genetic drift during branching leads to the mutation accumulation pattern that does not reflect the tree topology. Our modelling framework can help interpret and provide further insights into the empirical findings of genetic mosaicism in long-lived trees.

Molecular motors are responsible for carrying cellular transport of various membranous vesicles or organelles along cytoskeletal tracks. Transport of cellular cargos require high forces that are generated by motors working in groups. Hence, the properties of cargo transport can be modulated by varying various parameters such as cargo size and shape, microtubule geometry, motor number and their arrangement on cargo surface. Only those motors which are present in the contact zone on cargo surface have potential to bind to microtubule. Although earlier studies revealed the importance of cargo size, total motors attached to microtubule and their arrangement on cargo transport, yet how the contact zone influences binding of motors to microtubule largely remains unexplored. Here, it has been shown that contact zone is elliptical in shape for a spherical cargo and increases with cargo size for Kinesin-1 motors. To further understand the combined effect of elliptical contact zone and microtubule geometry on cargo transport, 3D mean-field model with uniform and clustered arrangement of motors for different cargo sizes and motor number has been used. Our findings indicate that cylindrical microtubule geometry maximizes the microtubule-bound motors which enhances the runlength and velocity of cargo transport. Our results show that microtubule-bound motors decrease with cargo size for uniform arrangement of motors on cargo thus decreasing its runlength and velocity, whereas in clustered arrangement, the number of microtubule-bound motors increase with cargo size which leads to increase in runlength and velocity.

The SARS-CoV-2 coronavirus continues to evolve with scores of mutations of the spike, membrane, envelope, and nucleocapsid structural proteins that impact pathogenesis. Infection data from nasal swabs, nasal PCR assays, upper respiratory samples, ex vivo cell cultures and nasal epithelial organoids reveal extreme variabilities in SARS-CoV-2 RNA titers within and between the variants. Some variabilities are naturally prone to clinical testing protocols and experimental controls. Here we focus on nasal viral load sensitivity arising from the timing of sample collection relative to onset of infection and from heterogeneity in the kinetics of cellular infection, uptake, replication, and shedding of viral RNA copies. The sources of between-variant variability are likely due to SARS-CoV-2 structural protein mutations, whereas within-variant population variability is likely due to heterogeneity in cellular response to that particular variant. With the physiologically faithful, agent-based mechanistic model of inhaled exposure and infection from (Chen et al., 2022), we perform statistical sensitivity analyses of the progression of nasal viral titers in the first 0-48 h post infection, focusing on three kinetic mechanisms. Model simulations reveal shorter latency times of infected cells (including cellular uptake, viral RNA replication, until the onset of viral RNA shedding) exponentially accelerate nasal viral load. Further, the rate of infectious RNA copies shed per day has a proportional influence on nasal viral load. Finally, there is a very weak, negative correlation of viral load with the probability of infection per virus-cell encounter, the model proxy for spike-receptor binding affinity.