We study an eco-evolutionary dynamics in finite populations of two haploid asexually reproducing allelic types. We focus on the quasi-neutral case when individual types differ only in their intrinsic birth and death rates but have the same expected lifetime reproductive output. We assume that the population size can fluctuate stochastically. We solve the Kolmogorov forward equation in the population whose size fluctuates only minimally and show that the fixation probability is decreasing with the increasing turnover rate. We also show that when the mutant's turnover is small enough, selection favors the mutant replacing residents. Similarly, when the turnover is high enough, selection opposes the replacement. This basic result has previously been demonstrated numerically for the contact process and shown analytically for the Moran process; the current paper extends this analysis to provide an analytical proof for the contact process. We also demonstrate numerically that our results extend for general fluctuating populations and beyond the quasi-neutral case.

Based on the physicochemical indexes of 20 amino acids and the Hungarian algorithm, each amino acid was mapped into a vector. And, the protein sequence can be represented as time series in eleven-dimensional space. In addition, the DTW algorithm was applied to calculate the distance between two time series to compare the similarities of protein sequences. The validity and accuracy of this method was illustrated by similarity comparison of ND5 proteins of nine species. Furthermore, homology analysis of eleven ACE2 proteins, which included human, Malayan pangolin and six species of bats, confirmed that the human had shorter evolutionary distance from the pangolin than those bats. The phylogenetic tree of spike protein sequences of 36 coronaviruses, which were divided into five groups, Class I, Class II, Class III, SARS-CoVs and COVID-19, was constructed.

Competition is a pervasive interaction known to structure ecological communities. The Lotka-Volterra (LV) model has been foundational for our understanding of competition, and trait-based LV models have been used to model community assembly and eco-evolutionary phenomena like diversification. The intrinsic growth rate function is determined by the underlying resource distribution and is a key determinant of the resulting diversity, traits and abundances of species. In these models, the width of the resource distribution relative to the width of the competition kernel has been identified as a key parameter that leads to diversification. However, studies have only investigated the impact of width at just a few discrete values, while also often assuming the intrinsic growth rate function to be unimodal. Thus, the impact of the underlying resource distribution's width and shape together remains incompletely explored, particularly for large, diverse communities. In this study, we vary its width continuously for two shapes (unimodal and bimodal) to explore its impact on community structure. When the resource distribution is very narrow in both the unimodal bimodal cases, competition is strong, leading to exclusion of all but the best-adapted species. Wider resource distributions allow stable coexistence, where the traits of the species depend on the shape of the resource distribution. Extremely wide resource distributions support a diverse community, where the strength of competition ultimately determines the diversity and traits of coexisting species, but their abundances reflect the underlying resource distribution. Further, competition acts to maximize the use of available resources among the competing species. For large communities, the shape of resource distribution becomes immaterial and the width determines the diversity. These results affirm and extend our understanding of limiting similarity.

The vascular endothelial growth factor (VEGF) receptor (VEGFR) system plays a role in cancer and many other diseases. It is widely accepted that VEGFR receptors dimerise in response to VEGF binding. However, analysis of these mechanisms and their implications for drug development still requires further exploration. In this paper, we present a mathematical model representing the binding of VEGF to VEGFR and the subsequent ligand-induced dimerisation. A key factor in this work is the qualitative and quantitative effect of binding cooperativity, which describes the effect that the binding of a ligand to a receptor has on the binding of that ligand to a second receptor, and the dimerisation of these receptors. We analyse the ordinary differential equation system at equilibrium, giving analytical solutions for the total amount of ligand bound. For time-course dynamics, we use numerical methods to explore possible behaviours under various parameter regimes, while perturbation analysis is used to understand the intricacies of these behaviours. Our simulation results show an excellent fit to experimental data, towards validating the model.

Theoretical studies have attempted to predict the effects of losing species and their relative importance to community-level processes, but little attention has been paid to specific interaction types such as omnivory and the trophic level of extinct species. Here, we use complex food web models to assess the importance of omnivores to community-level stability. We simulated food webs with varying size and complexity so we could remove omnivores from different trophic levels and track the fates of communities. Results show that food webs were more sensitive to the loss of omnivores from higher trophic levels, aside from the effects of network size and complexity on food web stability. These results underline the importance of omnivory by introducing dynamic switching across multiple energy paths. We expect this paper to shed light on additional factors that are fundamental when modelling extinctions such as the trophic level of lost species.

In group foraging animals, vigilance tends to decrease as group size increases. A forager in a group receives a vigilance benefit not only when it is being vigilant itself but also when a group mate is being vigilant. The many eyes hypothesis supposes that individuals exhibit lower vigilance in larger groups because of this. However, changes in safety resulting from the vigilance benefit conferred by group mates can change the decision to join or leave a group so as vigilance changes because of changes in group size, group size may also change in response to changes in vigilance. Additionally, individuals may have poor information about the vigilance strategies of their neighbors. We present a game theoretical model of vigilance that incorporates dynamic group sizes and does not require behavioral monitoring of the vigilance strategies of others. For systems at equilibrium, maximum vigilance decreases with increased group size. Furthermore, by varying intraspecific competition we show an inverse relationship between group size and vigilance. Thus, we provide a mechanism in support of the many eyes hypothesis from an evolutionary game theory perspective and conclude that variation in intraspecific competition and its effect on group size may be responsible for the relationship.

Intra-tumour heterogeneity (ITH) has a strong impact on the efficacy of the immune response against solid tumours. The number of sub-populations of cancer cells expressing different antigens and the percentage of immunogenic cells (i.e. tumour cells that are effectively targeted by immune cells) in a tumour are both expressions of ITH. Here, we present a spatially explicit stochastic individual-based model of the interaction dynamics between tumour cells and CD8+ T cells, which makes it possible to dissect out the specific impact of these two expressions of ITH on anti-tumour immune response. The set-up of numerical simulations of the model is defined so as to mimic scenarios considered in previous experimental studies. Moreover, the ability of the model to qualitatively reproduce experimental observations of successful and unsuccessful immune surveillance is demonstrated. First, the results of numerical simulations of this model indicate that the presence of a larger number of sub-populations of tumour cells that express different antigens is associated with a reduced ability of CD8+ T cells to mount an effective anti-tumour immune response. Secondly, the presence of a larger percentage of tumour cells that are not effectively targeted by CD8+ T cells may reduce the effectiveness of anti-tumour immunity. Ultimately, the mathematical model presented in this paper may provide a framework to help biologists and clinicians to better understand the mechanisms that are responsible for the emergence of different outcomes of immunotherapy.

The central role of metabolism in cell functioning and adaptation has given rise to countless studies on the evolution of enzyme-coding genes and network topology. However, very few studies have addressed the question of how enzyme concentrations change in response to positive selective pressure on the flux, considered a proxy of fitness. In particular, the way cellular constraints, such as resource limitations and co-regulation, affect the adaptive landscape of a pathway under selection has never been analyzed theoretically. To fill this gap, we developed a model of the evolution of enzyme concentrations that combines metabolic control theory and an adaptive dynamics approach, and integrates possible dependencies between enzyme concentrations. We determined the evolutionary equilibria of enzyme concentrations and their range of neutral variation, and showed that they differ with the properties of the enzymes, the constraints applied to the system and the initial enzyme concentrations. Simulations of long-term evolution confirmed all analytical and numerical predictions, even though we relaxed the simplifying assumptions used in the analytical treatment.

Insufficient testing capacity has been a critical bottleneck in the worldwide fight against COVID-19. Optimizing the deployment of limited testing resources has therefore emerged as a keystone problem in pandemic response planning. Here, we use a modified SEIR model to optimize testing strategies under a constraint of limited testing capacity. We define pre-symptomatic, asymptomatic, and symptomatic infected classes, and assume that positively tested individuals are immediately moved into quarantine. We further define two types of testing. Clinical testing focuses only on the symptomatic class. Non-clinical testing detects pre- and asymptomatic individuals from the general population, and a concentration parameter governs the degree to which such testing can be focused on high infection risk individuals. We then solve for the optimal mix of clinical and non-clinical testing as a function of both testing capacity and the concentration parameter. We find that purely clinical testing is optimal at very low testing capacities, supporting early guidance to ration tests for the sickest patients. Additionally, we find that a mix of clinical and non-clinical testing becomes optimal as testing capacity increases. At high but empirically observed testing capacities, a mix of clinical testing and non-clinical testing, even if extremely unfocused, becomes optimal. We further highlight the advantages of early implementation of testing programs, and of combining optimized testing with contact reduction interventions such as lockdowns, social distancing, and masking.

Seed formation is part of the reproductive cycle, leading to the accumulation of resistance stages that can withstand harsh environmental conditions for long periods of time. At the community level, multiple species with such long-lasting life stages can be more likely to coexist. While the implications of this process for biodiversity have been studied in terrestrial plants, seed banks are usually neglected in phytoplankton multispecies dynamic models, in spite of widespread empirical evidence for such seed banks. In this study, we build a metacommunity model of interacting phytoplankton species, including a resting stage supplying the seed bank. The model is parameterized with empirically-driven growth rate functions and field-based interaction estimates, which include both facilitative and competitive interactions. Exchanges between compartments (coastal pelagic cells, coastal resting cells on the seabed, and open ocean pelagic cells) are controlled by hydrodynamical parameters to which the sensitivity of the model is assessed. We consider two models, i.e., with and without a saturating effect of the interactions on the growth rates. Our results are consistent between models, and show that a seed bank allows to maintain all species in the community over 30 years. Indeed, a fraction of the species are vulnerable to extinction at specific times within the year, but this process is buffered by their survival in their resting stage. We thus highlight the potential role of the seed bank in the recurrent re-invasion of the coastal community, and of coastal environments in re-seeding oceanic regions. Moreover, the seed bank enables populations to tolerate stronger interactions within the community as well as more severe changes to the environment, such as those predicted in a climate change context. Our study therefore shows how resting stages may help phytoplanktonic diversity maintenance.

Computational models of biological processes provide one of the most powerful methods for a detailed analysis of the mechanisms that drive the behavior of complex systems. Logic-based modeling has enhanced our understanding and interpretation of those systems. Defining rules that determine how the output activity of biological entities is regulated by their respective inputs has proven to be challenging. Partly this is because of the inherent noise in data that allows multiple model parameterizations to fit the experimental observations, but some of it is also due to the fact that models become increasingly larger, making the use of automated tools to assemble the underlying rules indispensable. We present several Boolean function metrics that provide modelers with the appropriate framework to analyze the impact of a particular model parameterization. We demonstrate the link between a semantic characterization of a Boolean function and its consistency with the model's underlying regulatory structure. We further define the properties that outline such consistency and show that several of the Boolean functions under study violate them, questioning their biological plausibility and subsequent use. We also illustrate that regulatory functions can have major differences with regard to their asymptotic output behavior, with some of them being biased towards specific Boolean outcomes when others are dependent on the ratio between activating and inhibitory regulators. Application results show that in a specific signaling cancer network, the function bias can be used to guide the choice of logical operators for a model that matches data observations. Moreover, graph analysis indicates that commonly used Boolean functions become more biased with increasing numbers of regulators, supporting the idea that rule specification can effectively determine regulatory outcome despite the complex dynamics of biological networks.

Organisms consist of several genetic factors differing between species. However, the evolutionary effects of gene interactions on the evolutionary rate, adaptation, and divergence of organisms remain unknown. In a previous study, the 2-dimensional genetic algorithm (2DGA) program, including a gene interaction parameter, could simulate punctuated equilibrium under the disparity mode. Following this, we verified the effect of the number of gene interactions (gene cluster size) on evolution speed, adaptation, and divergence using the advanced 2DGA program. In this program, the population was replicated, mutated, and selected for 200,000 generations, and the fitness score, divergence, number of population, and genotype were output and plotted. The genotype data were used for evaluating the phylogenetic relations among the population. The gene cluster size 1) affected the disparity and parity mutagenesis modes differently, 2) determined the growth/exclusion rate and error threshold, and 3) accelerated or decelerated the population's speed of evolutionary advancement. In particular, when the gene cluster size expanded, the rate of increase in fitness scores decreased independently of the mutation rate and mode of mutation (disparity mode/parity mode). The mutation rate at the error threshold was also decreased by expanding the gene cluster size. Dendrograms traced the genotypes of the simulated population, indicating that the disparity mode caused the evolutionary process to enter 1) a stun mode, 2) an evolution mode, or 3) a divergence mode based on the mutation rate and gene cluster size, while the parity mode did not cause the population to enter a stun mode. Based on the above findings, we compared the predictions of the present study with evolution observed in the laboratory or the natural world and the processes of ongoing virus evolution, suggesting that our findings possibly explained the real evolution.

A general theory for competitive dynamics among many strains at the epidemiological level is required to understand polymorphisms in virulence, transmissibility, antibiotic resistance and other biological traits of infectious agents. Mathematical coinfection models have addressed specific systems, focusing on the criteria leading to stable coexistence or competitive exclusion, however, due to their complexity and nonlinearity, analytical solutions in coinfection models remain rare. Here we study a 2-strain Susceptible-Infected-Susceptible (SIS) compartmental model with co-infection/co-colonization, incorporating five strain fitness dimensions under the same framework: variation in transmissibility, duration of carriage, pairwise susceptibilities to coinfection, coinfection duration, and transmission priority effects from mixed coinfection. Taking advantage of a singular perturbation approach, under the assumption of strain similarity, we expose how strain dynamics on a slow timescale are explicitly governed by a replicator equation which encapsulates all traits and their interplay. This allows to predict explicitly not only the final epidemiological outcome of a given 2-player competition, but moreover, their entire frequency dynamics as a direct function of their relative variation and of strain-transcending global parameters. Based on mutual invasion fitnesses, we analyze and report rigorous results on transition phenomena in the 2-strain system, strongly mediated via endemic coinfection prevalence. We show that coinfection is not always a promoter of coexistence; instead, its effect to favour or prevent polymorphism is non-monotonic and depends on the type and level of phenotypic differentiation between strains. This framework offers a deeper analytical understanding of 2-strain competitive games in coinfection, with theoretical and practical applications in epidemiology, ecology and evolution.

The stress-gradient hypothesis (SGH) in ecology predicts that the strength and frequency of positive interspecific interactions, including processing chain commensalisms (PCCs), increase with environmental stress. Although observed in some empirical PCC studies, a recent theoretical study of PCCs using a consumer-resource-type model found that, given the model's assumptions, the SGH never occurs. To investigate if this is a true reflection of PCCs or merely an artefact of the model, in this study, we modified this earlier model formulation by incorporating generalized, monotonically increasing resource uptake functions in place of linear functions, and added a term to represent the spontaneous leakage of the downstream resource to the environment. Mathematical analyses of the model revealed two key insights: 1) the specific algebraic forms of the functional responses of the species in a PCC do not affect the long-term behaviour of the system; 2) the SGH can occur in a facilitative interaction only if the consumer-independent leakage rate of the downstream resource exceeds the consumer-independent input rate. The first insight shows that the outcomes of consumer-resource interactions are robust to details of the functional responses when the functional responses are monotonically increasing, while the second insight shows that the SGH is not a universal feature of positive interactions but instead holds only under a well-defined set of conditions which may vary between PCC interactions and the environmental contexts in which they take place.

In this paper, we conduct theoretical analyses on inferring the structure of gene regulatory networks. Depending on the experimental method and data type, the inference problem is classified into 20 different scenarios. For each scenario, we discuss the problem that with enough data, under what assumptions, what can be inferred about the structure. For scenarios that have been covered in the literature, we provide a brief review. For scenarios that have not been covered in literature, if the structure can be inferred, we propose new mathematical inference methods and evaluate them on simulated data. Otherwise, we prove that the structure cannot be inferred.

The Delta-Notch-Hes signaling pathway is involved in various developmental processes ranging from the formation of somites to the dynamic fine-grained patterns of cell types in developing or regenerating tissues. Such broad patterning capabilities rely in part on the versatile and tunable dynamics of the Notch-Hes feedback circuit eliciting both pulsatile and switching behaviors. This raises the theoretical issue of which specific spatiotemporal features emerge from lateral inhibition between cells that can display and transit between monostable, oscillatory and bistable regimes. To address this issue, I consider a discrete cell lattice model where intracellular dynamics is described by a phase-like variable and displays a typical cross-shaped phase diagram. Model analysis determines how the existence and stability of many spatially inhomogeneous and temporally synchronized states depends on key intracellular and intercellular parameters. It reveals a parameter-dependent multistability between those diverse spatiotemporal patterns, giving rise to tunable and robust developmental transition scenarios ensuring defect-free spatial patterns. Such broad repertoire and multistability of spatiotemporal patterns is corroborated with regulatory network modeling of the Delta-Notch-Hes pathway.

A recent study of adaptive dynamics of lysis propensity in temperate phages suggested that full lysogeny emerges as the outcome of bacteriophage evolution in a simple host-phage system. The conclusion is based on the premise that mutant strains necessarily appear in equilibrium host-phage environments. Revisiting the model, we show that the ecological system exhibits richer asymptotic dynamics and that, in a certain parameter regime, evolution may in fact drive lysis propensity towards an evolutionary singularity in which a non-zero proportion of phages initiate infection in a lytic cycle. These singularities act as points of evolutionary diversification, leading to periodic coexistence of two distinct phage strains on the evolutionary time-scale. One of the two strains in the dimorphic evolutionary singularity is fully lysogenic (in the sense that cell infection always leads to lysogeny), while the other is partially lytic. Our study thus highlights the importance of ecological interactions as a driver of evolution.

Temporal environmental stochasticity (TES), along with the variations of demographic rates associated with it, is ubiquitous in nature. Here we study the effect of TES on the species richness of diverse communities. In such communities the biodiversity at equilibrium reflects the balance between the rate at which new types are added (via migration, mutation or speciation) and the rate of extinction. We analyze a few generic models in which the speciation rate is fixed and TES affects the rate of extinction, and identify three different mechanisms. First, TES increases abundance variations and shortens extinction times, thus decreasing the species richness (destabilizing effect). Second, TES blurs the time-independent fitness differences between species, making the dynamics more symmetric and thereby increasing the diversity (neutralizing effect). Third, the storage effect allows TES to facilitate the invasion of inferior species, again contributing to the species richness. The stabilizing effect of storage declines significantly in diverse communities and it can overcome the destabilizing effect of TES only when environmental fluctuations are rapid enough.

This paper develops a detailed model of human and bovine erythrocytes quantifying the dependence of total cell volume upon composition of an aqueous solution in which it is immersed. The cytoplasm is represented as an aqueous solution of hemoglobin and salt (KCl or NaCl). Model-based analysis of literature data on human erythrocytes, and of new experiments with bovine erythrocytes, leads to two findings. First, the Boyle-van't Hoff plot for human erythrocytes is found to be well described based on a solid volume fraction of ∼0.3 in complete agreement with desiccation experiments. The linear portion of the calculated curve turns out to be numerically indistinguishable from the commonly used ideal model parameterized with an apparent osmotically inactive volume fraction of ∼0.5. This mathematical outcome explains the longstanding perceived (but actually nonexistent) disconnect between the aforementioned fractions ∼0.3 and ∼0.5. A corollarial implication is that the actual volume fraction of osmotically nonparticipant (vicinal) water is very small (∼0.035). Second, an initial crenation of bovine erythrocytes (which occurs in classical techniques for measuring membrane permeability) is found to increase their fragility to an extent which correlates well with the crenated cell volume, and would affect the permeability determination.

tBregs are a newly discovered subcategory of B regulatory cells, which are generated by breast cancer, resulting in the increase of Tregs and therefore in the death of NK cells. In this study, we use a mathematical and computational approach to investigate the complex interactions between the aforementioned cells as well as CD8+ T cells, CD4+ T cells and B cells. Furthermore, we use data fitting to prove that the functional response regarding the lysis of breast cancer cells by NK cells has a ratio-dependent form. Additionally, we include in our model the concentration of rituximab - a monoclonal antibody that has been suggested as a potential breast cancer therapy - and test its effect, when the standard, as well as experimental dosages, are administered.

This paper introduces a simple 1-dimensional map-based model of spiking neurons. During the past decades, dynamical models of neurons have been used to investigate the biology of human nervous systems. The models simulate experimental records of neurons' voltages using difference or differential equations. Difference neuronal models have some advantages besides the differential ones. They are usually simpler, and considering the cost of needed computations, they are more efficient. In this paper, a simple 1-dimensional map-based model of spiking neurons is introduced. Sample entropy is applied to analyze the complexity of the model's dynamics. The model can generate a wide range of time series with different firing rates and different levels of complexities. Besides, using some tools like bifurcation diagrams and cobwebs, the introduced model is analyzed.

Many models of within-host malaria infection dynamics have been formulated since the pioneering work of Anderson et al. in 1989. Biologically, the goal of these models is to understand what governs the severity of infections, the patterns of infectiousness, and the variation thereof across individual hosts. Mathematically, these models are based on dynamical systems, with standard approaches ranging from K-compartments ordinary differential equations (ODEs) to delay differential equations (DDEs), to capture the relatively constant duration of replication and bursting once a parasite infects a host red blood cell. Using malariatherapy data, which offers fine-scale resolution on the dynamics of infection across a number of individual hosts, we compare the fit and robustness of one of these standard approaches (K-compartments ODE) with a partial differential equations (PDEs) model, which explicitly tracks the "age" of an infected cell. While both models perform quite similarly in terms of goodness-of-fit for suitably chosen K, the K-compartments ODE model particularly overestimates parasite densities early on in infections when the number of repeated compartments is not large enough. Finally, the K-compartments ODE model (for suitably chosen K) and the PDE model highlight a strong qualitative connection between the density of transmissible parasite stages (i.e., gametocytes) and the density of host-damaging (and asexually-replicating) parasite stages. This finding provides a simple tool for predicting which hosts are most infectious to mosquitoes -vectors of Plasmodium parasites- which is a crucial component of global efforts to control and eliminate malaria.

We present a model of the hydraulic power required by the network of veins in a leaf with pinnate venation. The pinnate networks we study are dendritic networks with a single midrib and L levels of hierarchy, where L=2 corresponds to secondary veins branching from the midrib, L=3 additionally has tertiary veins branching from secondary veins, L=4 additionally has quaternary veins branching from tertiary veins, etc.We begin by utilizing the classic results of Murray to show that the minimal power required in a pipe of constant radius depends only on the length of the pipe and the volume flow rate to the 2/3 power. After then showing that the power required by the midrib is essentially independent of the number of secondary and higher order veins, we provide an explicit formula for the minimal total power required with L levels of hierarchy. The critical parameters in this model are the height and width of the leaf and the density ρ of vein terminations. We show how ρ can be estimated from published data on leaf vein density, and that for a very wide range of ρ the value of L minimizing the total power required is either 3 or 4. That is, three or four levels of leaf vein hierarchy suffice to minimize the required hydraulic power.

Trade is a complex, multi-faceted process that can contribute to the spread and persistence of disease. We here develop novel mechanistic models of supply. Our model is framed within a livestock trading system, where farms form and end trade partnerships with rates dependent on current demand, with these trade partnerships facilitating trade between partners. With these time-varying, stock dependent partnership and trade dynamics, our trading model goes beyond current state of the art modelling approaches. By studying instantaneous shocks to farm-level supply and demand we show that behavioural responses of farms lead to trading systems that are highly resistant to shocks with only temporary disturbances to trade observed. Individual adaptation in response to permanent alterations to trading propensities, such that animal flows are maintained, illustrates the ability for farms to find new avenues of trade, minimising disruptions imposed by such alterations to trade that common modelling approaches cannot adequately capture. In the context of endemic disease control, we show that these adaptations hinder the potential beneficial reductions in prevalence such changes to trading propensities have previously been shown to confer. Assessing the impact of a common disease control measure, post-movement batch testing, highlights the ability for our model to measure the stress on multiple components of trade imposed by such control measures and also highlights the temporary and, in some cases, the permanent disturbances to trade that post-movement testing has on the trading system.

Tuberculosis (TB), caused by infection with Mycobacterium tuberculosis (Mtb), is one of the world's deadliest infectious diseases and remains a significant global health burden. TB disease and pathology can present clinically across a spectrum of outcomes, ranging from total sterilization of infection to active disease. Much remains unknown about the biology that drives an individual towards various clinical outcomes as it is challenging to experimentally address specific mechanisms driving clinical outcomes. Furthermore, it is unknown whether numbers of immune cells in the blood accurately reflect ongoing events during infection within human lungs. Herein, we utilize a systems biology approach by developing a whole-host model of the immune response to Mtb across multiple physiologic and time scales. This model, called HostSim, tracks events at the cellular, granuloma, organ, and host scale and represents the first whole-host, multi-scale model of the immune response following Mtb infection. We show that this model can capture various aspects of human and non-human primate TB disease and predict that biomarkers in the blood may only faithfully represent events in the lung at early time points after infection. We posit that HostSim, as a first step toward personalized digital twins in TB research, offers a powerful computational tool that can be used in concert with experimental approaches to understand and predict events about various aspects of TB disease and therapeutics.

The microcirculation is the site of direct oxygen (O2) transfer from blood to tissue, and also of O2 delivery control via regulation of local blood flow. In addition, a number of diseases including type II diabetes mellitus (DMII) and sepsis are known to produce microcirculatory dysfunction in their early phases. Given the complexity of microvascular structure and physiology, and the difficulty of measuring tissue oxygenation at the micro-scale, mathematical modelling has been necessary for understanding the physiology and pathophysiology of O2 transport in the microcirculation and for interpreting in vivo experiments. To advance this area, a model of blood-tissue O2 transport in skeletal muscle was recently developed which uses continuously distributed capillaries and includes O2 diffusion, convection, and consumption. The present work extends this model to two adjacent layers of skeletal muscle with different blood flow rates and applies it to study steady-state O2 transport when flow regulation is stimulated using an O2 exchange chamber. To generate a model which may be validated through in vivo experiments, an overlying O2 permeable membrane is included. The model is solved using traditional methods including separation of variables and Fourier decomposition, and to ensure smooth profiles at the muscle-muscle and muscle-membrane interfaces, matching conditions are developed. The study presents qualitative verification for the model, using visualizations of tissue O2 pressure (PO2) distributions for varying capillary density (CD), and presents capillary velocity response values in the near layer for varying chamber PO2 under the assumption that outlet capillary O2 saturation (SO2) is equalized between adjacent layers. These compensatory velocity profiles, along with effective 'no-flux' chamber PO2 values, are presented for varying CD and tissue O2 consumption values. Insights gained from the two-layer model provide guidance for interpreting and planning future in vivo experiments, and also provide motivation for further development of the model to improve understanding of the interaction between O2 transport and blood flow regulation.