Update of bioRxiv. 2020 Nov 17;:

We discuss the dynamical robustness of biological networks represented by directed graphs, such as neural networks and gene regulatory networks. The theoretical results indicate that networks with low indegree variance and high outdegree variance are dynamically robust. We propose a machine learning method that gives equilibrium states to input-output networks with a recurrent hidden layer. We verify the theory by using the learned networks having various indegree and outdegree distributions. We also show that the basin of attraction of an equilibrium state is narrow when networks are dynamically robust.

Parasites can alter the behavior of animals. Such alterations could be a byproduct of infection or actively controlled and directed by the parasite. Ants infected with zombie ant fungi (Ophiocordyceps sp.) show behavioral changes culminating in the ant dying while biting into vegetation. To investigate the influence of the parasite on behavioral changes, we created an agent-based model that provides a prediction of how fungal infected ants move before death. The model shows how alterations in movement, such as an increased turning rate, within the normal range of ant behavior, can lead a host from the nest to the underside of a leaf. This demonstrates the simplicity in how such behavioral changes could evolve, as the fungal parasite could benefit from the natural behavior of the host, contesting a hypothesis of highly directed manipulation.

Social networks are clustered networks with short mean path length. In this work we analyze the disease dynamics in a class of this type of small-world networks composed of set of households and a set of workplaces. Individuals from each household are randomly assigned to workplaces. In both environments we assumed complete mixing and therefore we obtain highly clustered networks with short mean path lengths. Basic reproduction numbers were computed numerically and we show that at endemic equilibrium the average susceptible proportion <S/N> is different from the inverse of the basic reproduction number (R0-1). Therefore exist an exponent p≠1 for which <S/N>p=R0-1. Using this exponent we developed a mean field model which closely capture the disease dynamics in the network. Finally we outline how this model could be use to model vector-borne diseases in social networks.

Recent preclinical studies have shown that interstitial fluid pressure (IFP) within tumors can be heterogeneous Andersen et al. (2019). In that study tumors of two xenograft models, respectively, HL-16 cervical carcinoma and Panc-1 pancreatic carcinoma, were investigated. Significant heterogeneity in IFP was reported and it was proposed that this was associated with division of tissue into compartments separated by thick connective tissue bands for the HL-16 tumors and with dense collagen-rich extracellular matrix for the Panc-1 tumors. The purpose of the current work is to explore these experimental observations by using in silico generated tumor models. We consider a mathematical multiphase model which accounts for tumor cells, fibroblasts and interstitial fluid. The model has been trained to comply with experimental in vitro results reported in Shieh et al. (2011) which has identified autologous chemotaxis, ECM remodeling, and cell-fibroblast interaction as drivers for invasive tumor cell behavior. The in silico model is informed with parameters that characterize the leaky intratumoral vascular network, the peritumoral lymphatics which collect the fluid, and the density of ECM as represented through the hydraulic conductivity of the interstitial space. Heterogeneous distribution of solid stress may result in heterogeneous compression of blood vessels and, thus, heterogeneous vascular density inside the tumor. To mimic this we expose the in silico tumor to an intratumoral vasculature whose net effect of density of blood vesssels and vessel wall conductivity is varied through a 2D Gaussian variogram constrained such that the resulting IFPs lie within the range as reported from the preclinical study. The in silico cervical carcinoma model illustrates that sparse ECM was associated with uniform intratumoral IFP in spite of heterogeneous microvascular network, whereas compartment structures resulted in more heterogeneous IFP. Similarly, the in silico pancreatic model shows that heterogeneity in the microvascular network combined with dense ECM structure prevents IFP to even out and gives rise to heterogeneous IFP. The computer model illustrates how a heterogeneous invasive front might form where groups of tumor cells detach from the primary tumor and form isolated islands, a behavior which is natural to associate with metastatic propensity. However, unlike experimental studies, the current version of the in silico model does not show an association between metastatic propensity and elevated IFP.

Cancer stem cells (CSCs) often switch on their self-renewal programming aggressively to cause a relapse of cancer. Intriguingly, glucose triggers the proliferation propensities in CSCs by controlling the expression of the key transcription factor-like Nanog. However, the factors that critically govern this glucose-stimulated proliferation dynamics of CSCs remain elusive. Herein, by proposing a mathematical model of glucose-mediated Nanog regulation, we showed that the differential proliferation behavior of CSCs and cell-type similar to CSCs can be explained by considering the experimentally observed varied expression levels of key positive (STAT3) and negative (p53) regulators of Nanog. Our model reconciles various experimental observations and predicts ways to fine-tune the proliferation dynamics of these cell types in a context-dependent manner. In future, these modeling insights will be useful in developing improved therapeutic strategies to get rid of harmful CSCs.

This study develops a hemodynamic model involving the atrium, ventricle, veins, and arteries that can be calibrated to experimental results. It is a Windkessel model that incorporates an unsteady Bernoulli effect in the blood flow to the atrium. The model is represented by ordinary differential equations in terms of blood volumes in the compartments as state variables and it demonstrates the use of conductance instead of resistance to capture the effect of a non-leaking heart valve. The experimental results are blood volume data from 20 young (half of which are women) and 20 elderly (half of which are women) subjects during rest, inotropic stress (dobutamine), and chronotropic stress (glycopyrrolate). The model is calibrated to conform with data and physiological findings in 4 different levels. First, an optimization routine is devised to find model parameter values that give good fit between the model volume curves and blood volume data in the atrium and ventricle. Patient-specific information are used to get initial parameter values as a starting point of the optimization. Also, model pressure curves must show realistic behavior. Second, parametric bootstrapping is performed to establish the reliability of the optimal parameters. Third, statistical tests comparing mean optimal parameter values from young vs elderly subjects and women vs men are examined to support and present age and sex related differences in heart functions. Lastly, statistical tests comparing mean optimal parameter values from resting condition vs pharmacological stress are studied to verify and quantify the effects of dobutamine and glycopyrrolate to the cardiovascular system.

We study the influence of population heterogeneity on herd immunity level and on individual's vaccination decision making process. We first formulate the mathematical model in a population with two subgroups, based on different activity levels or different susceptibilities. The herd immunity threshold is derived and discussed. It is calculated that the required vaccine coverage level for herd immunity in a heterogeneous mixing population can be varied significantly. The required vaccine coverage level is lower than the classical herd immunity level, if the vaccine coverage level in the more active group or more susceptible group is higher than the other subgroup. It is suggested that the classical herd immunity levels can be misleading in the process of planning mass vaccination programs. The analysis is further extended to study the population with more subgroups. We then study the formal vaccination games to simulate the process of vaccination decision making, in either homogeneous or heterogeneous mixing populations. It is proved that the Nash equilibrium in the vaccination game is not unique if population heterogeneity is considered. Moreover, herd immunity is not achieved if individuals are solely driven by self-interests.

Intermediate host species provide a crucial link in the emergence of zoonotic infectious diseases, serving as a population where an emerging pathogen can mutate to become human-transmissible. Identifying such species is thus a key component of predicting and possibly mitigating future epidemics. Despite this importance, intermediate host species have not been investigated in much detail, and have generally only been identified by testing for the presence of pathogens in multiple candidate species. In this paper, we present a mathematical model able to identify likely intermediate host species for emerging zoonoses based on ecological data for the candidates and epidemiological data for the pathogen. Since coronaviruses frequently emerge through intermediate host species and, at the time of writing, pose an urgent pandemic threat, we apply the model to the three emerging coronaviruses of the twenty-first century, accurately predicting palm civets as intermediate hosts for SARS-CoV-1 and dromedary camels as intermediate hosts for MERS. Further, we suggest mink, pangolins, and ferrets as intermediate host species for SARS-CoV-2. With the capacity to evaluate intermediate host likelihood among different species, researchers can focus testing for possible infection sources and interventions more effectively.

African swine fever virus (ASFV) leads to a highly contagious, lethal and economically devastating disease among pigs. Since no effective treatment for the disease, it is crucial to investigate its transmission mechanism and control strategies in large-scale pig farms. We first established a toy model to explore ASFV spread in one pig unit. Then a switching patch model was developed to capture its spread from one initial epidemic pig house consecutively to others, even the whole farm. Assessing innocent culling rates of three large-scale epidemic pig farms in Jiangsu Province showed that it is unnecessary to slaughter all pigs in the farms compulsively. Then we explored how the disinfection and fixation of employees impact ASFV spread in the farms. To control ASFV, we can block or slow down its spreading by improving the efficiency of disinfection and decreasing employee population to some extend. We can also shrink potential areas to be infected by properly improving the matching refinement degree among employees and houses. Some essential requirements for large-scale pig farms are presented to reduce their ASFV spreading risk, which can be helpful for animal health authorities in establishing regulation to standardize large-scale pig farms.

Social animals often share information about the location of resources, such as a food source or a new nest-site. One well-studied communication strategy in ants is tandem running, whereby a leader guides a recruit to a resource. Tandem running is considered an example of animal teaching because a leader adjusts her behaviour and invests time to help another ant to learn the location of a resource more efficiently. Tandem running also has costs, such as waiting inside the nest for a leader and a reduced walking speed. Whether and when these costs outweigh the benefits of tandem running is not well understood. We developed an agent-based simulation model to investigate the conditions that favour communication by tandem running during foraging. We predicted that the spatio-temporal distribution of food sources, colony size and the ratio of scouts and recruits affect colony foraging success. Our results suggest that tandem running is favoured when food sources are hard to find, differ in energetic value and are long lasting. These results mirror the findings of simulations of honeybee communication. Scouts locate food sources faster than tandem followers in some environments, suggesting that tandem running may fulfil the criteria of teaching only in some situations. Furthermore, tandem running was only beneficial above a critical colony size threshold. Taken together, our model suggests that there is a considerable parameter range that favours colonies that do not use communication by tandem running, which could explain why many ants with small colony sizes forage solitarily.

Modularity is a common feature of the nervous system across species and scales. Although it has been qualitatively investigated in network science, very little is known about how it affects spike signal transmission in neuronal networks at the mesoscopic level. Here, a neuronal network model is built to simulate dynamic interactions among different modules of neuronal networks. This neuronal network model follows the organizational principle of modular structure. The neurons can generate spikes like biological neurons, and changes in the strength of synaptic connections conform to the STDP learning rule. Based on this neuronal network model, we first quantitatively studied whether and to what extent the connectivity within and between modules can affect spike signal transmission, and found that spike signal transmission heavily depends on the connectivity between modules, but has little to do with the connectivity within modules. More importantly, we further found that the spike activity of a module can be predicted according to the spike activities of its adjacent modules through building a resting-state functional connectivity matrix.

Growing genetic lineage mapping experiments have definitively shown a wide-ranging plasticity of vascular smooth muscle cells (VSMCs) in atherosclerotic plaque and suggested that VSMCs can modulate their phenotypes in response to plaque microenvironment. Here, a multiscale hybrid discrete-continuous (HDC) modeling system is established to investigate the complex role of VSMC phenotypic switching within atherosclerotic lesions. The cellular behaviors of VSMCs and macrophages, including proliferation, migration, phenotypic transformation and necrosis, are determined by cellular automata (CA) rules in discrete model. While the dynamics of plaque microenvironmental factors, such as lipid, extracellular matrix (ECM) and chemokines, are described by continuous reaction-diffusion equations in macroscopy. The simulation results demonstrate how the VSMC activities change the extracellular microenvironment and consequently affect the plaque morphology and stability. The regulation of VSMC phenotypes can affect not only the plaque morphology (necrotic core size and fibrous cap thickness) but also the deposition and distribution of microenvironmental factors (lipoprotein, ECM, and chemokines). In addition, it is found that plaque vulnerability can be inhibited by blocking VSMC transdifferentiation to a macrophage-like state and promoting it to a myofibroblastic phenotype, which suggests that targeting VSMC phenotypic switching could be a potential and promising therapeutic strategy for atherosclerosis.

In this study, we develop a methodology for model reduction and selection informed by global sensitivity analysis (GSA) methods. We apply these techniques to a control model that takes systolic blood pressure and thoracic tissue pressure data as inputs and predicts heart rate in response to the Valsalva maneuver (VM). The study compares four GSA methods based on Sobol' indices (SIs) quantifying the parameter influence on the difference between the model output and the heart rate data. The GSA methods include standard scalar SIs determining the average parameter influence over the time interval studied and three time-varying methods analyzing how parameter influence changes over time. The time-varying methods include a new technique, termed limited-memory SIs, predicting parameter influence using a moving window approach. Using the limited-memory SIs, we perform model reduction and selection to analyze the necessity of modeling both the aortic and carotid baroreceptor regions in response to the VM. We compare the original model to systematically reduced models including (i) the aortic and carotid regions, (ii) the aortic region only, and (iii) the carotid region only. Model selection is done quantitatively using the Akaike and Bayesian Information Criteria and qualitatively by comparing the neurological predictions. Results show that it is necessary to incorporate both the aortic and carotid regions to model the VM.

Elements within biological systems interact and frequently self-organize from initially disordered states into highly structured patterns. The local self-activation and lateral inhibition mechanism, derived from the coupling between two reacting and diffusing chemicals, has been believed to be one of the main causes for biological pattern formation. Graded positional information can be produced by the limited diffusion of one single signaling molecule through cell populations with no pre-patterns being required. We demonstrate, using multiscale computations, that spontaneous symmetry breaking can be driven within expanding and non-expanding cell populations, without local self-enhancement of activators and long-range inhibition. Instead, cells can self-organize into structured gene patterns via a combination of timing gene expression in cells and the graded positional information which has been coupled to the gene expression. We show that the genetic symmetry breaking in expanding E. coli populations occurs at a critical colony size, which is independent of the cell doubling time but scales with the diffusion speed of the signaling molecule. We also show the quasi-3D structure of gene patterns, and observe that the wave length of periodic genetic stripes is in proportion to the genetic oscillation cycle time and in inverse proportion to cell doubling time. Our results provide insights into relevant biological development processes.

In order to interpret animal behavior we need to understand how they see the world. As directly testing color discrimination in animals is difficult and time consuming, it is important to develop theoretical models based in the properties of visual systems. One of the most successful for the prediction of color discrimination behavior is the receptor noise-limited (RNL) model, which depends only on the level of noise in photoreceptors and opponent mechanisms. Here a complementary approach to model construction is used, and optimal color discrimination properties are obtained using information theoretical tools, for the early stages of visual systems. It is shown here that, for most biologically relevant conditions the optimal discrimination function of an ideal observer coincides with the one obtained with the RNL model. Furthermore, within this framework the influence of opponency can be studied by considering models with and without that mechanism but with exactly the same parameters at the level of photoreceptors. As an example, it is shown here that opponency is necessary to explain the discrimination of monochromatic stimuli in honeybees, but not in budgerigars. Since this is a consequence of the narrowing of absorption spectra of photoreceptors, produced by the presence of oil droplets, this could also be true for most other species of birds. This suggests that in order to study opponency in birds, stimuli should have a relatively wide spectrum.

Most traits expressed by organisms, such as gene expression profiles, developmental trajectories, behavioural sequences and reaction norms are function-valued traits (colloquially "phenotypically plastic traits"), since they vary across an individual's age and in response to various internal and/or external factors (state variables). Furthermore, most organisms live in populations subject to limited genetic mixing and are thus likely to interact with their relatives. We here formalise selection on genetically determined function-valued traits of individuals interacting in a group-structured population, by deriving the marginal version of Hamilton's rule for function-valued traits. This rule simultaneously gives a condition for the invasion of an initially rare mutant function-valued trait and its ultimate fixation in the population (invasion thus implies substitution). Hamilton's rule thus underlies the gradual evolution of function-valued traits and gives rise to necessary first-order conditions for their uninvadability (evolutionary stability). We develop a novel analysis using optimal control theory and differential game theory, to simultaneously characterise and compare the first-order conditions of (i) open-loop traits - functions of time (or age) only, and (ii) closed-loop (state-feedback) traits - functions of both time and state variables. We show that closed-loop traits can be represented as the simpler open-loop traits when individuals do not interact or when they interact with clonal relatives. Our analysis delineates the role of state-dependence and interdependence between individuals for trait evolution, which has implications to both life-history theory and social evolution.

Two recent high profile studies have attempted to use edge (branch) length ratios from large sets of phylogenetic trees to determine the relative ages of genes of different origins in the evolution of eukaryotic cells. This approach can be straightforwardly justified if substitution rates are constant over the tree for a given protein. However, such strict molecular clock assumptions are not expected to hold on the billion-year timescale. Here we propose an alternative set of conditions under which comparisons of edge length distributions from multiple sets of phylogenies of proteins with different origins can be validly used to discern the order of their origins. We also point out scenarios where these conditions are not expected to hold and caution is warranted.

During the outbreak of emerging infectious diseases, information dissemination dynamics significantly affects the individuals' psychological and behavioral changes, and consequently influences on the disease transmission. To investigate the interaction of disease transmission and information dissemination dynamics, we proposed a multi-scale model which explicitly models both the disease transmission with saturated recovery rate and information transmission to evaluate the effect of information transmission on dynamic behaviors. Considering time variation between information dissemination, epidemiological and demographic processes, we obtained a slow-fast system by reasonably introducing a sufficiently small quantity. We carefully examined the dynamics of proposed system, including existence and stability of possible equilibria and existence of backward bifurcation, by using the fast-slow theory and directly investigating the full system. We then compared the dynamics of the proposed system and the essential thresholds based on two methods, and obtained the similarity between the basic dynamical behaviors of the slow system and that of the full system. Finally, we parameterized the proposed model on the basis of the COVID-19 case data in mainland China and data related to news items, and estimated the basic reproduction number to be 3.25. Numerical analysis suggested that information transmission about COVID-19 pandemic caused by media coverage can reduce the peak size, which mitigates the transmission dynamics during the early stage of the COVID-19 pandemic.

Hepatitis B virus (HBV) with its high prevalence and death toll is one of the most important infectious diseases to study. Yet, there is very little progress in the development of within-host models for HBV, which has subsequently hindered our understanding of this virus. The uncertainty around the proliferation of infected hepatocytes has been studied but never in association with other important biological continuous events such as integrations and superinfections. This is despite the fact that these processes affect the diversity and composition of infected cell population in the liver and an improved understanding of the cellular composition will undoubtedly assist in strategizing against this viral infection. Here, we developed novel mathematical models that incorporate these key biological processes and analyzed them both analytically and numerically. Unaffected by the extent of integrated DNA (IDNA), the outcome of HBV infection was primarily dictated by the balance between processes generating and killing infected hepatocytes containing covalent closed circular DNA (cccDNA). The superinfection was found to be a key process in the spread of HBV infection as its exclusion could not reproduce experimentally observed composition of infected hepatocytes at peak of acute HBV infection, a stage where our model predicts that infected hepatocytes most likely carry both cccDNA and IDNA. Our analysis further suggested the existence of some form of selective advantage of infected hepatocytes containing only IDNA to explain the viral dynamics observed during antiviral treatment and the transition from peak to acute infection. Finally, the fine line between liver destruction and resolution of acute HBV infection was found to be highly influenced by the fate of cccDNA during cellular proliferation.

The genetic code consists in a set of rules used by living organisms to translate genomic information, contained in genes, into proteins; every amino acid is coded by a set of nucleotide triplets or codons. We refer to codon choice as the choice of a given codon, among the synonymous available ones, to code a given amino acid occurrence. The aim of this work is to shed light on the pivotal role that codon choice plays in regulating the timing of translation process, through patterns of low and high translation efficiency codons. A translation efficiency value, namely codon score, was associated to each codon through a formula based on the number of tRNAs gene copies able to translate the given codon. By using codon scores, those k-mers of the proteome of Saccharomyces cerevisiae, showing low and high average scores associated to the correspondent codons, were computed. The analysis of distribution of both low and high average score k-mers clearly showed that, in particular for higher k-mer size, they occur much more than expected, strongly suggesting a functional role. Moreover performed analysis highlighted that significant k-mers preferentially occur in some protein folding classes, such as those containing alpha helices, and in some functional classes mainly involved in transcription process while codon choice seems to have a very low impact in proteins associated to energy production and metabolism. The relationship between secondary structures and significant k-mers was investigated, revealing that low score k-mers tend to preferentially occur in coil or close to coil regions and almost never in beta sheets, while high score k-mers preferentially occur in alpha helices, avoiding beta sheets, and close to coil regions for high k-mer sizes. Finally the analysis of distribution of significant codon patterns along the proteins highlighted a relevant enrichment of low average score k-mers at the 5' end of protein-coding sequences in the region from 5th to 25th amino acid.

The paper presents an attempt to integrate the classical evolutionary game theory based on replicator dynamics and the state-based approach of Houston and McNamara. In the new approach, individuals have different heritable strategies; however, individuals carrying the same strategy can differ in terms of state, role or the situation in which they act. Thus, the classical replicator dynamics is completed by the additional subsystem of differential equations describing the dynamics of transitions between different states. In effect, the interactions described by game structure, in addition to the demographic payoffs (constituted by births and deaths), can lead to the change in state of the competing individuals. Special cases of reversible and irreversible incremental stage-structured models, where the state changes can describeenergy accumulation, developmental steps or aging, are derived for discrete and continuous versions. The new approach is illustrated using the example of the Owner-Intruder game with explicit dynamics of the role changes. The new model presents a generalization of the demographic version of the Hawk-Dove game,with the difference being that the opponents in the game are drawn from two separate subpopulations consisting of Owners and Intruders. Here, the Intruders check random nest sites and play the Hawk-Dove game with the Owner if they are occupied. Meanwhile, the Owners produce newborns that become Intruders, since they must find a free nest site to reproduce. An interesting feedback mechanism is produced via the fluxes of individuals between the different subpopulations. In addition, the population growth suppression mechanism resulting from the fixation Bourgeois strategy is analyzed.

Erratum in J Theor Biol. 2022 Feb 7;534:110959.

Repeated games are useful models to analyze long term interactions of living species and complex social phenomena. Zero-determinant (ZD) strategies in repeated games discovered by Press and Dyson in 2012 enforce a linear payoff relationship between a focal player and the opponent. This linear relationship can be set arbitrarily by a ZD player. Hence, a subclass of ZD strategies can fix the opponent's expected payoff and another subclass of the strategies can exceed the opponent for the expected payoff. Since this discovery, theories for ZD strategies are extended to cope with various natural situations. It is especially important to consider the theory of ZD strategies for repeated games with a discount factor and observation errors because it allows the theory to be applicable in the real world. Recent studies revealed their existence of ZD strategies even in repeated games with both factors. However, the conditions for the existence has not been sufficiently analyzed. Here, we mathematically analyzed the conditions in repeated games with both factors. First, we derived the thresholds of a discount factor and observation errors which ensure the existence of Equalizer and positively correlated ZD (pcZD) strategies, which are well-known subclasses of ZD strategies. We found that ZD strategies exist only when a discount factor remains high as the error rates increase. Next, we derived the conditions for the expected payoff of the opponent enforced by Equalizer as well as the conditions for the slope and base line payoff of linear lines enforced by pcZD. As a result, we found that, as error rates increase or a discount factor decreases, the conditions for the linear line that Equalizer or pcZD can enforce become strict.

Continuous Glucose Monitoring (CGM) produces long time-series of noisy observations of a single variable (tissue glucose concentration), whose evolution may be explained by a dynamical model. In order to represent the unknown mixture of possible control mechanisms of different orders affecting the measured variable, a fractional differential approach seems justified. In any case, variations in food intake and/or physical activity ought to be taken into account if a plausible interpretation of the dynamics is to be obtained. In the present work, the mathematical construction and the numerical implementation of a Fractional Differential Equations (FDE) initial value problem are systematically reviewed, with the intent of offering the reader a concise and mathematically rigorous description of this approach. An FDE model for CGM is formulated: the model includes compartments for stomach and intestinal glucose contents and for blood and tissue (subcutaneous) glucose concentrations, as well as the shock effects of food ingestion and of increased glucose consumption due to physical activity. The model parameters, including the (non-integer) order of differentiation, are estimated from CGM observations on six Type 1 diabetic patients. The best-fit fractional orders for the six subjects range from 1.59 to 2.13. For comparison, best fits have also been computed for all subjects using an average fractional order of 1.9 and integer orders of 1 and 2.The results indicate that in the case of CGM the fractional differential model, which should be physiologically more appropriate, in fact fits the data much better than the first-order model and also better than the 2nd-order model.

According to the competitive-exclusion principle, the number n of regulating variables describing a given community dynamics is an upper bound on the number of species (or types or morphs) that can coexist at equilibrium. On occasion, it is possible to reformulate a model with a lower number of regulating variables than appeared in the initial specification. We call the smallest number of such variables the dimension of the environmental feedback, or environmental dimension for short. For studying which species can invade a community, it is enough to know the sign of each species' long-term growth rate, i.e., invasion fitness. Therefore, different indicators of population growth - so-called fitness proxies, such as the basic reproduction number-are sometimes preferred. However, as we show, different fitness proxies may have different dimensions. Fundamental characteristics such as the environmental dimension should not depend on such arbitrary choices. Here, we resolve this difficulty by introducing a refined definition of environmental dimension that focuses on neutral fitness contours. On this basis, we show that this definition of environmental dimension is not only unambiguous, i.e., independent of the choice of fitness proxy, but also constructive, i.e., applicable without needing to assess an infinite number of possible fitness proxies. We then investigate how to determine environmental dimensions by analysing the two components of the environmental feedback: the impact map describing how a community's resident species affect the regulating variables and the sensitivity map describing how population growth depends on the regulating variables. The dimension of the impact map is lower than n when the set of feasible environments is of lower dimension than n, and the dimension of the sensitivity map is lower than n when not all n regulating variables affect the sign of population growth independently. While the minimum of the dimensions of the impact and sensitivity maps provides an upper bound on the environmental dimension, the combined effect of the two maps can result in an even lower environmental dimension, which happens when the sensitivity map is insensitive to some aspects of the impact map's image. To facilitate the applications of the framework introduced here, we illustrate all key concepts with detailed worked examples. In view of these results, we claim that the environmental dimension is the ultimate generalization of the traditional and widely used notions of the "number of regulating variables" or the "number of limiting factors", and is thus the sharpest generally applicable upper bound on the number of species that can robustly coexist in a community.

Although mortality increases with age in most organisms, senescence is missing from models of parasite evolution. Since virulence evolves according to the host's mortality, and since virulence influences the intensity of transmission, which determines the average age at infection and thus the mortality rate of a senescing host, we expected that epi-evolutionary feedbacks would underlie the evolution of virulence in a population of senescing hosts. We tested this idea by extending an age-structured model of epidemiological dynamics with the parasite's evolution. A straightforward prediction of our model is that stronger senescence forces the evolution of higher virulence. However, the model also reveals that the evolved virulence depends on the average age at infection, giving an evolutionary feedback with the epidemiological situation, a prediction not found when assuming a constant mortality rate with age. Additionally, and in contrast to most models of parasite evolution, we found that the virulence at the evolutionary equilibrium is influenced by whether the force of infection depends on the density or on the frequency of infected hosts, due to changes in the average age at infection. Our findings suggest that ignoring age-specific effects, and in particular senescence, can give misleading predictions about parasite evolution.

In a population of interacting individuals, the environment for interactions often changes due to individuals' behaviors, which in turn drive the evolution of individuals' behaviors. The interplay between the environment and individuals' behaviors has been demonstrated to remarkably influence the evolutionary outcomes. In reality, in highly cognitive species such as social primates and human beings, individuals are often capable of perceiving the environment change and then differentiate their strategies across different environment states. We propose a model of environmental feedback with state-dependent strategies: individuals have perceptions of distinct environment states and therefore take distinct sub-strategies under each of them; based on the sub-strategy, individuals then decide their behaviors; their behaviors subsequently modify the environment state. We use the theory of stochastic games and evolutionary dynamics to analyze this idea. We find that when environment changes slower than behaviors, state-dependent strategies (i.e. taking different sub-strategies under different environment states) can outperform state-independent strategies (i.e. taking an identical sub-strategy under all environment states), such as Win-Stay, Lose-Shift, the most leading strategy among state-independent strategies. The intuition is that delayed environmental feedback provides chances for individuals with state-dependent strategies to exploit those with state-independent strategies. Our results hold (1) in both well-mixed and structured populations; (2) when the environment switches between two or more states. Furthermore, the environment changing rate decides if state-dependent strategies benefit global cooperation. The evolution sees the rise of the cooperation level for fast environment switching and the decrease otherwise. Our work stresses that individuals' perceptions of different environment states are beneficial to their survival and social prosperity in a changing world.

To be able to deal with uncertainty is of primary importance to most living organisms. When cues provide information about the state of the environment, organisms can use them to respond flexibly. Life forms have evolved complex adaptations and sensory mechanisms to use these environmental cues and extract valuable information about the environment. Previous work has shown a theoretical limit to the amount of fitness benefit possible to be extracted from the cues. We show that the previously used information theoretical approaches can be generalised to scenarios involving any potential relationship between the number of possible phenotypes and environmental states. Such cases are relevant when physiological constraints or complex ecological scenarios lead to the number of environmental states exceeding potential phenotypes. We illustrate cases in which these scenarios can emerge: along environmental gradients, such as geographical transects or complex environments, where organisms adopt different bet-hedging strategies, switching stochastically between phenotypes or developing intermediate ones. In conclusion, we develop an information-theoretic extensible approach for investigating and quantifying fitness in ecological studies.

Circadian clocks allow living organisms to anticipate and adapt to the daily variations of the environment. The interlocked feedback loops of the transcription factors network in the plant clock generate oscillations with expression peaks at specific times of the day. In this work, we explore the effect of molecular noise on the behavior of the plant circadian clock through numerical simulations. The influence of system size, photoperiod, and mutations of clock genes on the robustness of the oscillations are discussed. Our simulations show that the oscillations remain robust when the mRNA and protein levels are in the range of a few hundreds molecules. Entrainment by light-dark cycles enhances the robustness compared to constant conditions. Multiple light inputs and inter-cellular coupling also contribute to the robustness of the oscillations. The comparison between deterministic and stochastic simulations of single and double mutants shows that stochasticity does not qualitatively affect the behaviour of mutants but that they do not have the same robustness to noise. Finally, the model shows that noise can induce transitions between two limit cycles in a birhythmic clock mutant.

The maintenance of the proliferative cell niche is critical to epithelial tissue morphology and function. In this paper we investigate how current modelling methods can result in the erroneous loss of proliferative cells from the proliferative cell niche. Using an established model of the inter-follicular epidermis we find there is a limit to the proliferative cell densities that can be maintained in the basal layer (the niche) if we do not include additional mechanisms to stop the loss of proliferative cells from the niche. We suggest a new methodology that enables maintenance of a desired homeostatic population of proliferative cells in the niche: a rotational force is applied to the two daughter cells during the mitotic phase of division to enforce a particular division direction. We demonstrate that this new methodology achieves this goal. This methodology reflects the regulation of the orientation of cell division.

Capelin (Mallotus villosus) is a short-lived (1-4 years) fish species, that plays a crucial role by dominating the intermediate trophic level in the Barents Sea. Several episodes of extreme biomass decline (collapse) have been observed during the last three decades. We postulate that these collapses might be regulated by food availability (bottom-up effect) and/or by time discrepancy between capelin feeding and abundance of its prey (match-mismatch hypothesis). This paper investigates our postulate using a model consisting of a set of coupled differential equations to describe the predator-prey system, with a single delay term, τ, in description of the predator dynamics. We derive theoretical conditions on τ, as well as determine how changes in these conditions define different stability regimes of the system. Unconstrained optimization is used to calculate optimal model parameters by fitting the predator-prey model to empirical data. The optimization results are combined with those from the theoretical analysis, to make inference about the empirical system stability. Our results show that Hopf bifurcation occurs in the predatory-prey system when τ exceeds a theoretically derived value τ∗>0. This value represents the critical time for prey availability in advance of the optimal predator growth period.Set into an ecological context, our findings provide mathematical evidence for validity of the match-mismatch hypothesis and a bottom-up effect for capelin.

Phagocytosis is a complex process by which phagocytes such as lymphocytes or macrophages engulf and destroy foreign bodies called pathogens in a tissue. The process is triggered by the detection of antibodies that trigger signaling mechanisms that control the changes of the cellular cytoskeleton needed for engulfment of the pathogen. A mathematical model of the entire process would be extremely complicated, because the signaling and cytoskeletal changes produce large mechanical deformations of the cell. Recent experiments have used a confinement technique that leads to a process called frustrated phagocytosis, in which the membrane does not deform, but rather, signaling triggers actin waves that propagate along the boundary of the cell. This eliminates the large-scale deformations and facilitates modeling of the wave dynamics. Herein we develop a model of the actin dynamics observed in frustrated phagocytosis and show that it can replicate the experimental observations. We identify the key components that control the actin waves and make a number of experimentally-testable predictions. In particular, we predict that diffusion coefficients of membrane-bound species must be larger behind the wavefront to replicate the internal structure of the waves. Our model is a first step toward a more complete model of phagocytosis, and provides insights into circular dorsal ruffles as well.

Learning is thought to be achieved by the selective, activity dependent, adjustment of synaptic connections. Individual learning can also be very hard and/or slow. Social, supervised, learning from others might amplify individual, possibly mainly unsupervised, learning by individuals, and might underlie the development and evolution of culture. We studied a minimal neural network model of the interaction of individual, unsupervised, and social supervised learning by communicating "agents". Individual agents attempted to learn to track a hidden fluctuating "source", which, linearly mixed with other masking fluctuations, generated observable input vectors. In this model data are generated linearly, facilitating mathematical analysis. Learning was driven either solely by direct observation of input data (unsupervised, Hebbian) or, in addition, by observation of another agent's output (supervised, Delta rule). To make learning more difficult, and to enhance biological realism, the learning rules were made slightly connection-inspecific, so that incorrect individual learning sometimes occurs. We found that social interaction can foster both correct and incorrect learning. Useful social learning therefore presumably involves additional factors some of which we outline.

Prions are proteins that cause fatal neurodegenerative diseases. The misfolded conformation adopted by prions can be transmitted to other normally folded proteins. Therapeutics to stop prion proliferation have been studied experimentally; however, it is not clear how the combination of different types of treatments can decrease the growth rate of prions in the brain. In this article, we combine the implementation of pharmacological chaperones and interferons to develop a novel model using a non-linear system of ordinary differential equations and study the quantitative effects of these two treatments on the growth rate of prions. This study aims to identify how the two treatments affect prion proliferation, both individually and in tandem. We analyze the model, and qualitative global results on the disease-free and disease equilibria are proved analytically. Numerical simulations, using parameter values from in vivo experiments that provide a pharmaceutically important demonstration of the effects of these two treatments, are presented here. This mathematical model can be used to identify and optimize the best combination of the treatments within their safe ranges.

The effects of partner choice have been documented in a large number of biological systems such as sexual markets, interspecific mutualisms, or human cooperation. There are, however, a number of situations in which one would expect this mechanism to play a role, but where no such effect has ever been demonstrated. This is the case in particular in many intraspecific interactions, such as collective hunts, in non-human animals. Here we use individual-based simulations to solve this apparent paradox. We show that the conditions for partner choice to operate are in fact restrictive. They entail that individuals can compare social opportunities and choose the best. The challenge is that social opportunities are often rare because they necessitate the co-occurrence of (i) at least one available partner, and (ii) a resource to exploit together with this partner. This has three consequences. First, partner choice cannot lead to the evolution of cooperation when resources are scarce, which explains that this mechanism could never be observed in many cases of intraspecific cooperation in animals. Second, partner choice can operate when partners constitute in themselves a resource, which is the case in sexual interactions and interspecific mutualisms. Third, partner choice can lead to the evolution of cooperation when individuals live in a rich environment, and/or when they are highly efficient at extracting resources from their environment.

We study the influence of spatial heterogeneity on the antiviral activity of mouse embryonic fibroblasts (MEF) infected with influenza A. MEF of type Ube1L-/- are composed of two distinct sub-populations, the strong type that sustains a strong viral infection and the weak type, sustaining a weak viral load. We present new data on the virus load infection of Ube1L-/-, which have been micro-printed in a checker board pattern of different sizes of the inner squares. Surprisingly, the total viral load at one day after inoculation significantly depends on the sizes of the inner squares. We explain this observation by using a reaction diffusion model and we show that mathematical homogenization can explain the observed inhomogeneities. If the individual patches are large, then the growth rate and the carrying capacity will be the arithmetic means of the patches. For finer and finer patches the average growth rate is still the arithmetic mean, however, the carrying capacity uses the harmonic mean. While fitting the PDE to the experimental data, we also predict that a discrepancy in virus load would be unobservable after only half a day. Furthermore, we predict the viral load in different inner squares that had not been measured in our experiment and the travelling distance the virions can reach after one day.

Mound structures that soil termites build have diverse morphologies. Previous observational studies documented that mounds are built to provide regulated environments for the termites that live within them and their structures are formed in ways to support this purpose under the influence of the mounds' immediate environment. The objective of this study is to provide a methodology and a predictive computational model to investigate the reason behind the different but systematic shapes of termite mounds, considering all the relevant forces imposed on them and their thermoregulatory and gas-exchange functions. The gas-exchange function accounts for the capacity of the mound to diffuse metabolic gases generated in the mound's underground nest, while the thermoregulatory function satisfies the connection between the underground nest and deep ground temperatures. The proposed predictive model is based on the principles of heat transfer and thermodynamics and allows optimized mechanically stable structures to freely emerge. The results indicate that, while the model is free to generate any mechanically stable structure, under the relevant environmental and metabolic conditions, it produces structures with forms and geometrical characteristics similar to those of natural mounds. Investigation of the connection between the local environment and the mound shapes indicated that the Sun and wind play an important role in the mound structural form. Mounds exposed to stronger solar irradiance exhibit cone-shaped structures that are pointed towards the Sun, while shaded mounds are observed to be vertical domes. The local wind is observed to affect the external shape of the mound by preventing them to grow tall while controlling the features of the internal structure. By investigating the similarities between structures in different regions (i.e., India, Namibia, and Brazil), it is revealed that, unlike mounds with a strong need for gas-exchange, mounds with a significant demand for thermoregulation exhibit deeper nests, thicker external walls, and well-defined cone- (as opposed to the dome-) shaped structures. Overall, the form of termite mounds is strongly correlated to both regulatory functions and local environments, and the resulting mound shape arises as a combination of these factors.

The epidemiology of sexually transmitted infections (STIs) is inherently linked to host mating dynamics. Studies across many taxa show that adult sex ratio, a major determinant of host mating dynamics, is often skewed - sometimes strongly - toward males or females. However, few predictions exist for the effects of skewed sex ratio on STI epidemiology, and none when coupled with sex biased disease characteristics. Here we use mathematical modelling to examine how interactions between sex ratio and disease characteristics affect STI prevalence in males and females. Notably, we find that while overall disease prevalence peaks at equal sex ratios, prevalence per sex peaks at skewed sex ratios. Furthermore, disease characteristics, sex-biased or not, drive predictable differences in male and female STI prevalence as sex ratio varies, with higher transmission and lower virulence generally increasing differences between the sexes for a given sex ratio. Our work reveals new insights into how STI prevalence in males and females depends on a complex interaction between host population sex ratio and disease characteristics.

Bovine Viral Diarrhea (BVD) is a cattle disease that causes substantial financial losses, in particular to the dairy industry. Hence, several countries including Germany introduced compulsory disease control programs. For the case of Germany in particular, all animals had to be tested and persistently infected animals (PI animals) were removed from the population. The program was successful in reducing the number of PI animals, but was overtly expensive. Alternative approaches were therefore discussed to eliminate the remaining PI animals and alter the testing system in order to reduce costs. Contributing to these efforts, we developed an agent-based model that aimed to cover all relevant aspects of the disease biology and would allow to evaluate different control strategies. For the biological part of the infection spread, the model includes horizontal and vertical transmission, transient and persistent infections. Moreover, several control strategies including import of animals, trade restrictions, vaccination, as well as various testing schemes were included. The model was furthermore defined to be stochastic, event-driven and hierarchical, with cattle movements as the main route of spreading between farms. For the spread within farms, we included susceptible-infected-recovered (SIR) dynamics with an additional permanently infectious class. The interaction between the farms was described by a supply and demand farm manager mechanism governing the network structure and dynamics. Additionally, we carried out a sensitivity analysis of the input parameters to study the impact of extreme values on the model. Since the population size in the model is limited, we tested the influence of the initial population size on the model results. Our results showed that the model could accurately describe the dynamics of the disease in the presence and absence of disease control. Although we developed the model for the spread of BVD, it may be adapted to similar diseases of cattle and swine.

It is likely that RNA replication began non-enzymatically, and that polymerases were later selected to speed up the process. We consider replication mechanisms in modern viruses and ask which of these is possible non-enzymatically, using mathematical models and experimental data found in the literature to estimate rates of RNA synthesis and replication. Replication via alternating plus and minus strands is found in some single-stranded RNA viruses. However, if this occurred non-enzymatically it would lead to double-stranded RNA that would not separate. With some form of environmental cycling, such as temperature, salinity, or pH cycling, double-stranded RNA can be melted to form single-stranded RNA, although re-annealing of existing strands would then occur much faster than synthesis of new strands. We show that re-annealing blocks this form of replication at a very low concentration of strands. Other kinds of viruses synthesize linear double strands from single strands and then make new single strands from double strands via strand-displacement. This does not require environmental cycling and is not blocked by re-annealing. However, under non-enzymatic conditions, if strand-displacement occurs from a linear template, we expect the incomplete new strand to be almost always displaced by the tail end of the old strand through toehold-mediated displacement. A third kind of replication in viruses and viroids is rolling-circle replication which occurs via strand-displacement on a circular template. Rolling-circle replication does not require environmental cycling and is not prevented by toehold-mediated displacement. Rolling-circle replication is therefore expected to occur non-enzymatically and is a likely starting point for the evolution of polymerase-catalysed replication.

Development of resistance to chemotherapy in cancer patients strongly effects the outcome of the treatment. Due to chemotherapeutic agents, resistance can emerge by Darwinian evolution. Besides this, acquired drug resistance may arise via changes in gene expression. A recent discovery in cancer research uncovered a third possibility, indicating that this phenotype conversion can occur through the transfer of microvesicles from resistant to sensitive cells, a mechanism resembling the spread of an infectious agent. We present a model describing the evolution of sensitive and resistant tumour cells considering Darwinian selection, Lamarckian induction and microvesicle transfer. We identify three threshold parameters which determine the existence and stability of the three possible equilibria. Using a simple Dulac function, we give a complete description of the dynamics of the model depending on the three threshold parameters. We also establish an agent based model as a spatial version of the ODE model and compare the outputs of the two models. We find that although the ODE model does not provide spatial information about the structure of the tumour, it is capable to determine the outcome in terms of tumour size and distribution of cell types. We demonstrate the possible effects of increasing drug concentration, and characterize the possible bifurcation sequences. Our results show that the presence of microvesicle transfer cannot ruin a therapy that otherwise leads to extinction, however it may doom a partially successful therapy to failure.

Sea turtles complete migrations across vast distances, covering entire ocean basins. To track these migrations, satellite tracking tags are attached to their shells. The impact of these tags must be considered to ensure that turtles' natural behavior is not artificially and adversely impacted through tag-related drag, and that the data collected by a small sample of sea turtles accurately represents the larger population. Additionally, it can be difficult to study animal energetics in the field over large migration distances. In this work, we modify a computational behavior model to study how satellite tracking tags affect turtle migration behavior. Our agent based model contains synthetic magnetic field environments that are used for navigation cues, an ocean current, resource distributions that represent locations of food, and an agent that attempts to migrate to several different goals. The agent loses energy as it progresses, and searches for the resource distributions to replenish itself. Our novel simulation framework demonstrates the relationship between an agent's available energy capacity, its energy consumption based on mechanical power expended, and its ability to navigate to all migratory goal points. This study can be utilized to (1) probe the impacts of an animal's energy capacity and foraging behavior on its resulting navigation and ecology, (2) guide future satellite tag designs, and (3) develop usage recommendations for a suitable tracking tag based on the type of experiment being conducted. Our model can be expanded beyond sea turtles to study other marine species (e.g., sharks, whales). Additionally, this model could be expanded to other domains within the marine environment. For example, it could be modified to examine design trade-offs in remotely operated vehicles (ROVs), which share many of the same operational constraints as sea turtles and other migratory species.

Maintaining cerebral blood flow is critical for adequate neuronal function. Previous computational models of brain capillary networks have predicted that heterogeneous cerebral capillary flow patterns result in lower brain tissue partial oxygen pressures PO2). However, these previous models have often considered simple capillary networks in terms of their geometric properties. In this current work, we developed and analyzed computational models of brain capillary networks to determine how perturbations of network properties impact tissue oxygen levels. The models include variabilities in both their geometric (segment lengths and diameters) and three-dimensional, topological structure. Two classes of capillary network models are considered. The first consists of equations for the oxygen partial pressure, PO2, in both a capillary network and the surrounding tissue. In order to gain insight into the behavior of this detailed model, we also consider a reduced model for changes in PO2 in just the capillary network. The main result is that for a general class of networks, random perturbations of either segment diameters or conductances will always, on average, decrease the average tissue oxygen levels. This result is supported through both simulations of the models and mathematical analysis. Our results promise to expand our understanding of cerebral capillary blood flow and its impact on the brain function in health and disease.

In cancer, treatment failure and disease recurrence have been associated with small subpopulations of cancer cells with a stem-like phenotype. In this paper, we develop and investigate a phenotype-structured model of solid tumour growth in which cells are structured by a stemness level, which varies continuously between stem-like and terminally differentiated behaviours. Cell evolution is driven by proliferation and death, as well as advection and diffusion with respect to the stemness structure variable. Here, the magnitude and sign of the advective flux are allowed to vary with the oxygen level. We use the model to investigate how the environment, in particular oxygen levels, affects the tumour's population dynamics and composition, and its response to radiotherapy. We use a combination of numerical and analytical techniques to quantify how under physiological oxygen levels the cells evolve to a differentiated phenotype and under low oxygen level (i.e., hypoxia) they de-differentiate. Under normoxia, the proportion of cancer stem cells is typically negligible and the tumour may ultimately become extinct whereas under hypoxia cancer stem cells comprise a dominant proportion of the tumour volume, enhancing radio-resistance and favouring the tumour's long-term survival. We then investigate how such phenotypic heterogeneity impacts the tumour's response to treatment with radiotherapy under normoxia and hypoxia. Of particular interest is establishing how the presence of radio-resistant cancer stem cells can facilitate a tumour's regrowth following radiotherapy. We also use the model to show how radiation-induced changes in tumour oxygen levels can give rise to complex re-growth dynamics. For example, transient periods of hypoxia induced by damage to tumour blood vessels may rescue the cancer cell population from extinction and drive secondary regrowth.

The quality of hosts for a parasitoid wasp may be influenced by attributes such as host size or species, with high quality for successful development usually coincident with high quality for larger offspring. This is not always the case: for the Scelionid wasp Trissolcus basalis, oviposition in eggs of the Brown Marmorated Stink Bug, Halyomorpha halys, rather than of the normal host, the Southern Green Stink Bug, Nezara viridula, leads to lower offspring survival, but survivors can be unusually large. Adult female T. basalis engage in contests for host access. As larger contestants are typically favoured in contests between parasitoids, the larger size of surviving offspring may compensate for the mortality of others. We construct a general game-theoretic model to explore whether size advantage can sustain a maternal preference to utilize a more deadly host species. We find that size advantage alone is unlikely to sustain a shift in host preference, yet such an outcome is possible when size asymmetries act simultaneously with advantages in host possession (ownership effect). Halyomorpha halys is an invasive pest of major agro-economic importance in Europe and the Americas, and use of its eggs as hosts by native parasitoids such as T. basalis has been seen as an evolutionary trap due to their high developmental mortality. Our model suggests that the recently discovered effect of host choice on offspring size may provide an escape from the trap via effects on contest biology of T. basalis which could foster a more stable association with H. halys. An evolutionary shift in the reproductive value of H. halys could increase the efficiency of T. basalis as a biological control agent of this invasive stink bug pest.

In this work, we study the behavior of a time-delayed mutually repressive auto-activating three-gene system. Delays are introduced to account for the location difference between DNA transcription that leads to production of messenger RNA and its translation that result in protein synthesis. We study the dynamics of the system using numerical simulations, computational bifurcation analysis and mathematical analysis. We find Hopf bifurcations leading to stable and unstable rotation in the system, and we study the rotational behavior as a function of cyclic mutual repression parameter asymmetry between each gene pair in the network. We focus on how rotation co-exists with a stable heteroclinic flow linking the three saddles in the system. We find that this coexistence allows for a transition between two markedly different types of rotation leading to strikingly different phenotypes. One type of rotation belongs to Hopf-induced rotation while the other type, belongs to heteroclinic cycling between three saddle nodes in the system. We discuss the evolutionary and biological implications of our findings.