Hybrid models of chemotaxis with application to leukocyte migration.

Many chemical and biological systems involve reacting species with vastly different numbers of molecules/agents. Hybrid simulations model such phenomena by combining discrete (e.g., agent-based) and continuous (e.g., partial differential equation- or PDE-based) descriptors of the dynamics of reactants with small and large numbers of molecules/agents, respectively. We present a stochastic hybrid algorithm to model a stage of the immune response to inflammation, during which leukocytes reach a pathogen via chemotaxis. While large numbers of chemoattractant molecules justify the use of a PDE-based model to describe the spatiotemporal evolution of its concentration, relatively small numbers of leukocytes and bacteria involved in the process undermine the veracity of their continuum treatment by masking the effects of stochasticity and have to be treated discretely. Motility and interactions between leukocytes and bacteria are modeled via random walk and a stochastic simulation algorithm, respectively. Since the latter assumes the reacting species to be well mixed, the discrete component of our hybrid algorithm deploys stochastic operator splitting, in which the sequence of the diffusion and reaction operations is determined autonomously during each simulation step. We conduct a series of numerical experiments to ascertain the accuracy and computational efficiency of our hybrid simulations and, then, to demonstrate the importance of randomness for predicting leukocyte migration and fate during the immune response to inflammation.

Efficient models of polymerization applied to FtsZ ring assembly in Escherichia coli.

High protein concentrations complicate modeling of polymer assembly kinetics by introducing structural complexity and a large variety of protein forms. We present a modeling approach that achieves orders of magnitude speed-up by replacing distributions of lengths and widths with their average counterparts and by introducing a hierarchical classification of species and reactions into sets. We have used this model to study FtsZ ring assembly in Escherichia coli The model's prediction of key features of the ring formation, such as time to reach the steady state, total concentration of FtsZ species in the ring, total concentration of monomers, and average dimensions of filaments and bundles, are all in agreement with the experimentally observed values. Besides validating our model against the in vivo observations, this study fills some knowledge gaps by proposing a specific structure of the ring, describing the influence of the total concentration in short and long kinetics processes, determining some characteristic mechanisms in polymer assembly regulation, and providing insights about the role of ZapA proteins, critical components for both positioning and stability of the ring.

Quantification of Predictive Uncertainty in Models of FtsZ ring assembly in Escherichia coli.

Quantitative predictions of FtsZ protein polymerization are essential for understanding the self-regulating mechanisms in biochemical systems. Due to structural complexity and parametric uncertainty, existing kinetic models remain incomplete and their predictions error-prone. To address such challenges, we perform probabilistic uncertainty quantification and global sensitivity analysis of the concentrations of various protein species predicted with a recent FtsZ protein polymerization model. Our results yield a ranked list of modeling shortcomings that can be improved in order to develop more accurate predictions and more realistic representations of key mechanisms of such biochemical systems and their response to changes in internal or external conditions. Our conclusions and improvement recommendations can be extended to other kinetics models.

Effects of Hydraulic Soil Properties on Vegetation Pattern Formation in Sloping Landscapes.

Current models of vegetation pattern formation rely on a system of weakly nonlinear reaction-diffusion equations that are coupled by their source terms. While these equations, which are used to describe a spatiotemporal planar evolution of biomass and soil water, qualitatively capture the emergence of various types of vegetation patterns in arid environments, they are phenomenological and have a limited predictive power. We ameliorate these limitations by deriving the vertically averaged Richards' equation to describe flow (as opposed to "diffusion") of water in partially saturated soils. This establishes conditions under which this nonlinear equation reduces to its weakly nonlinear reaction-diffusion counterpart used in the previous models, thus relating their unphysical parameters (e.g., diffusion coefficient) to the measurable soil properties (e.g., hydraulic conductivity) used to parameterize the Richards equation. Our model is valid for both flat and sloping landscapes and can handle arbitrary topography and boundary conditions. The result is a model that relates the environmental conditions (e.g., precipitation rate, runoff and soil properties) to formation of multiple patterns observed in nature (such as stripes, labyrinth and spots).

Posttransfusion Increase of Hematocrit per se Does Not Improve Circulatory Oxygen Delivery due to Increased Blood Viscosity.

BACKGROUND: Blood transfusion is used to treat acute anemia with the goal of increasing blood oxygen-carrying capacity as determined by hematocrit (Hct) and oxygen delivery (DO2). However, increasing Hct also increases blood viscosity, which may thus lower DO2 if the arterial circulation is a rigid hydraulic system as the resistance to blood flow will increase. The net effect of transfusion on DO2 in this system can be analyzed by using the relationship between Hct and systemic blood viscosity of circulating blood at the posttransfusion Hct to calculate DO2 and comparing this value with pretransfusion DO2. We hypothesized that increasing Hct would increase DO2 and tested our hypothesis by mathematically modeling DO2 in the circulation. METHODS: Calculations were made assuming a normal cardiac output (5 L/min) with degrees of anemia ranging from 5% to 80% Hct deficit. We analyzed the effects of transfusing 0.5 or more units of 300 cc of packed red blood cells (PRBCs) at an Hct of 65% and calculated microcirculatory DO2 after accounting for increased blood viscosity and assuming no change in blood pressure. Our model accounts for O2 diffusion out of the circulation before blood arriving to the nutritional circulation and for changes in blood flow velocity. The immediate posttransfusion DO2 was also compared with DO2 after the transient increase in volume due to transfusion has subsided. RESULTS: Blood transfusion of up to 3 units of PRBCs increased DO2 when Hct (or hemoglobin) was 60% lower than normal, but did not increase DO2 when administered before this threshold. CONCLUSIONS: After accounting for the effect of increasing blood viscosity on blood flow owing to increasing Hct, we found in a mathematical simulation of DO2 that transfusion of up to 3 units of PRBCs does not increase DO2, unless anemia is the result of an Hct deficit greater than 60%. Observations that transfusions occasionally result in clinical improvement suggest that other mechanisms possibly related to increased blood viscosity may compensate for the absence of increase in DO2.

Shear-Induced Nitric Oxide Production by Endothelial Cells.

We present a biochemical model of the wall shear stress-induced activation of endothelial nitric oxide synthase (eNOS) in an endothelial cell. The model includes three key mechanotransducers: mechanosensing ion channels, integrins, and G protein-coupled receptors. The reaction cascade consists of two interconnected parts. The first is rapid activation of calcium, which results in formation of calcium-calmodulin complexes, followed by recruitment of eNOS from caveolae. The second is phosphorylation of eNOS by protein kinases PKC and AKT. The model also includes a negative feedback loop due to inhibition of calcium influx into the cell by cyclic guanosine monophosphate (cGMP). In this feedback, increased nitric oxide (NO) levels cause an increase in cGMP levels, so that cGMP inhibition of calcium influx can limit NO production. The model was used to predict the dynamics of NO production by an endothelial cell subjected to a step increase of wall shear stress from zero to a finite physiologically relevant value. Among several experimentally observed features, the model predicts a highly nonlinear, biphasic transient behavior of eNOS activation and NO production: a rapid initial activation due to the very rapid influx of calcium into the cytosol (occurring within 1-5 min) is followed by a sustained period of activation due to protein kinases.

Efficient Multiscale Models of Polymer Assembly.

Protein polymerization and bundling play a central role in cell physiology. Predictive modeling of these processes remains an open challenge, especially when the proteins involved become large and their concentrations high. We present an effective kinetics model of filament formation, bundling, and depolymerization after GTP hydrolysis, which involves a relatively small number of species and reactions, and remains robust over a wide range of concentrations and timescales. We apply this general model to study assembly of FtsZ protein, a basic element in the division process of prokaryotic cells such as Escherichia coli, Bacillus subtilis, or Caulobacter crescentus. This analysis demonstrates that our model outperforms its counterparts in terms of both accuracy and computational efficiency. Because our model comprises only 17 ordinary differential equations, its computational cost is orders-of-magnitude smaller than the current alternatives consisting of up to 1000 ordinary differential equations. It also provides, to our knowledge, a new insight into the characteristics and functioning of FtsZ proteins at high concentrations. The simplicity and versatility of our model render it a powerful computational tool, which can be used either as a standalone descriptor of other biopolymers' assembly or as a component in more complete kinetic models.

A boundary-layer solution for flow at the soil-root interface.

Transpiration, a process by which plants extract water from soil and transmit it to the atmosphere, is a vital (yet least quantified) component of the hydrological cycle. We propose a root-scale model of water uptake, which is based on first principles, i.e. employs the generally accepted Richards equation to describe water flow in partially saturated porous media (both in a root and the ambient soil) and makes no assumptions about the kinematic structure of flow in a root-soil continuum. Using the Gardner (exponential) constitutive relation to represent the relative hydraulic conductivities in the Richards equations and treating the root as a cylinder, we use a matched asymptotic expansion technique to derive approximate solutions for transpiration rate and the size of a plant capture zone. These solutions are valid for roots whose size is larger than the macroscopic capillary length of a host soil. For given hydraulic properties, the perturbation parameter used in our analysis relates a root's size to the macroscopic capillary length of the ambient soil. This parameter determines the width of a boundary layer surrounding the soil-root interface, within which flow is strictly horizontal (perpendicular to the root). Our analysis provides a theoretical justification for the standard root-scale cylindrical flow model of plant transpiration that imposes a number of kinematic constraints on water flow in a root-soil continuum.

Identifying transport behavior of single-molecule trajectories.

Models of biological diffusion-reaction systems require accurate classification of the underlying diffusive dynamics (e.g., Fickian, subdiffusive, or superdiffusive). We use a renormalization group operator to identify the anomalous (non-Fickian) diffusion behavior from a short trajectory of a single molecule. The method provides quantitative information about the underlying stochastic process, including its anomalous scaling exponent. The classification algorithm is first validated on simulated trajectories of known scaling. Then it is applied to experimental trajectories of microspheres diffusing in cytoplasm, revealing heterogeneous diffusive dynamics. The simplicity and robustness of this classification algorithm makes it an effective tool for analysis of rare stochastic events that occur in complex biological systems.

Hematocrit dispersion in asymmetrically bifurcating vascular networks.

Quantitative modeling of physiological processes in vasculatures requires an accurate representation of network topology, including vessel branching. We propose a new approach for reconstruction of vascular network, which determines how vessel bifurcations distribute red blood cells (RBC) in the microcirculation. Our method follows the foundational premise of Murray's law in postulating the existence of functional optimality of such networks. It accounts for the non-Newtonian behavior of blood by allowing the apparent blood viscosity to vary with discharge hematocrit and vessel radius. The optimality criterion adopted in our approach is the physiological cost of supplying oxygen to the tissue surrounding a blood vessel. Bifurcation asymmetry is expressed in terms of the amount of oxygen consumption associated with the respective tissue volumes being supplied by each daughter vessel. The vascular networks constructed with our approach capture a number of physiological characteristics observed in in vivo studies. These include the nonuniformity of wall shear stress in the microcirculation, the significant increase in pressure gradients in the terminal sections of the network, the nonuniformity of both the hematocrit partitioning at vessel bifurcations and hematocrit across the capillary bed, and the linear relationship between the RBC flux fraction and the blood flow fraction at bifurcations.

Non-Newtonian flow of blood in arterioles: consequences for wall shear stress measurements.

OBJECTIVE: Our primary goal is to investigate the effects of non-Newtonian blood properties on wall shear stress in microvessels. The secondary goal is to derive a correction factor for the Poiseuille-law-based indirect measurements of wall shear stress. METHODS: The flow is assumed to exhibit two distinct, immiscible and homogeneous fluid layers: an inner region densely packed with RBCs, and an outer cell-free layer whose thickness depends on discharge hematocrit. The cell-free layer is assumed to be Newtonian, while rheology of the RBC-rich core is modeled using the Quemada constitutive law. RESULTS: Our model provides a realistic description of experimentally observed blood velocity profiles, tube hematocrit, core hematocrit, and apparent viscosity over a wide range of vessel radii and discharge hematocrits. CONCLUSIONS: Our analysis reveals the importance of incorporating this complex blood rheology into estimates of WSS in microvessels. The latter is accomplished by specifying a correction factor, which accounts for the deviation of blood flow from the Poiseuille law.

Vegetation pattern formation due to interactions between water availability and toxicity in plant-soil feedback.

Development of a comprehensive theory of the formation of vegetation patterns is still in progress. A prevailing view is to treat water availability as the main causal factor for the emergence of vegetation patterns. While successful in capturing the occurrence of multiple vegetation patterns in arid and semiarid regions, this hypothesis fails to explain the presence of vegetation patterns in humid environments. We explore the rich structure of a toxicity-mediated model of the vegetation pattern formation. This model consists of three PDEs accounting for a dynamic balance between biomass, water, and toxic compounds. Different (ecologically feasible) regions of the model's parameter space give rise to stable spatial vegetation patterns in Turing and non-Turing regimes. Strong negative feedback gives rise to dynamic spatial patterns that continuously move in space while retaining their stable topology.

Anomalous diffusion of single particles in cytoplasm.

The crowded intracellular environment poses a formidable challenge to experimental and theoretical analyses of intracellular transport mechanisms. Our measurements of single-particle trajectories in cytoplasm and their random-walk interpretations elucidate two of these mechanisms: molecular diffusion in crowded environments and cytoskeletal transport along microtubules. We employed acousto-optic deflector microscopy to map out the three-dimensional trajectories of microspheres migrating in the cytosolic fraction of a cellular extract. Classical Brownian motion (BM), continuous time random walk, and fractional BM were alternatively used to represent these trajectories. The comparison of the experimental and numerical data demonstrates that cytoskeletal transport along microtubules and diffusion in the cytosolic fraction exhibit anomalous (nonFickian) behavior and posses statistically distinct signatures. Among the three random-walk models used, continuous time random walk provides the best representation of diffusion, whereas microtubular transport is accurately modeled with fractional BM.

Probability density function method for Langevin equations with colored noise.

Understanding the mesoscopic behavior of dynamical systems described by Langevin equations with colored noise is a fundamental challenge in a variety of fields. We propose a new approach to derive closed-form equations for joint and marginal probability density functions of state variables. This approach is based on a so-called large-eddy-diffusivity closure and can be used to model a wide class of non-Markovian processes described by the noise with an arbitrary correlation function. We demonstrate the accuracy of the proposed probability density function method for several linear and nonlinear Langevin equations.

Hypertonic treatment of acute respiratory distress syndrome.

Many viral infections, including the COVID-19 infection, are associated with the hindrance of blood oxygenation due to the accumulation of fluid, inflammatory cells, and cell debris in the lung alveoli. This condition is similar to Acute Respiratory Distress Syndrome (ARDS). Mechanical positive-pressure ventilation is often used to treat this condition, even though it might collapse pulmonary capillaries, trapping red blood cells and lowering the lung's functional capillary density. We posit that the hyperosmotic-hyperoncotic infusion should be explored as a supportive treatment for ARDS. As a first step in verifying the feasibility of this ARDS treatment, we model the dynamics of alveolar fluid extraction by osmotic effects. These are induced by increasing blood plasma osmotic pressure in response to the increase of blood NaCl concentration. Our analysis of fluid drainage from a plasma-filled pulmonary alveolus, in response to the intravenous infusion of 100 ml of 1.28 molar NaCl solution, shows that alveoli empty of fluid in approximately 15 min. These modeling results are in accordance with available experimental and clinical data; no new data were collected. They are used to calculate the temporal change of blood oxygenation, as oxygen diffusion hindrance decreases upon absorption of the alveolar fluid into the pulmonary circulation. Our study suggests the extraordinary speed with which beneficial effects of the proposed ARDS treatment are obtained and highlight its practicality, cost-efficiency, and avoidance of side effects of mechanical origin.

Autoregulation and mechanotransduction control the arteriolar response to small changes in hematocrit.

Here, we present an analytic model of arteriolar mechanics that accounts for key autoregulation mechanisms, including the myogenic response and the vasodilatory effects of nitric oxide (NO) in the vasculature. It couples the fluid mechanics of blood flow in arterioles with solid mechanics of the vessel wall and includes the effects of wall shear stress- and stretch-induced endothelial NO production. The model can be used to describe the regulation of blood flow and NO transport under small changes in hematocrit and to analyze the regulatory response of arterioles to small changes in hematocrit. Our analysis revealed that the experimentally observed paradoxical increase in cardiac output with small increases in hematocrit results from the combination of increased NO production and the effects of a strong myogenic response modulated by elevated levels of WSS. Our findings support the hypothesis that vascular resistance varies inversely with blood viscosity for small changes in hematocrit in a healthy circulation that responds to shear stress stimuli. They also suggest beneficial effects independent of changes in O(2) carrying capacity associated with the postsurgical transfusion of one or two units of blood.

PEG-albumin supraplasma expansion is due to increased vessel wall shear stress induced by blood viscosity shear thinning.

We studied the extreme hemodilution to a hematocrit of 11% induced by three plasma expanders: polyethylene glycol (PEG)-conjugated albumin (PEG-Alb), 6% 70-kDa dextran, and 6% 500-kDa dextran. The experimental component of our study relied on microelectrodes and cardiac output to measure both the rheological properties of plasma-expander blood mixtures and nitric oxide (NO) bioavailability in vessel walls. The modeling component consisted of an analysis of the distribution of wall shear stress (WSS) in the microvessels. Our experiments demonstrated that plasma expansion with PEG-Alb caused a state of supraperfusion with cardiac output 40% above baseline, significantly increased NO vessel wall bioavailability, and lowered peripheral vascular resistance. We attributed this behavior to the shear thinning nature of blood and PEG-Alb mixtures. To substantiate this hypothesis, we developed a mathematical model of non-Newtonian blood flow in a vessel. Our model used the Quemada rheological constitutive relationship to express blood viscosity in terms of both hematocrit and shear rate. The model revealed that the net effect of the hemodilution induced by relatively low-viscosity shear thinning PEG-Alb plasma expanders is to reduce overall blood viscosity and to increase the WSS, thus intensifying endothelial NO production. These changes act synergistically, significantly increasing cardiac output and perfusion due to lowered overall peripheral vascular resistance.

Impact of endothelium roughness on blood flow.

Cell free layer (CFL), a plasma layer bounded by the red blood cell (RBC) core and the endothelium, plays an important physiological role. Its width affects the effective blood viscosity as well as the scavenging and production of nitric oxide (NO). Measurements of the CFL and its spatio-temporal variability are highly uncertain, exhibiting random fluctuations. Yet traditional models of blood flow and NO scavenging treat the CFL's bounding surfaces as deterministic and smooth. We investigate the effects of the endothelium roughness and uncertain (random) spatial variability on blood flow and the estimates of effective blood viscosity.

Stochastic operator-splitting method for reaction-diffusion systems.

Many biochemical processes at the sub-cellular level involve a small number of molecules. The local numbers of these molecules vary in space and time, and exhibit random fluctuations that can only be captured with stochastic simulations. We present a novel stochastic operator-splitting algorithm to model such reaction-diffusion phenomena. The reaction and diffusion steps employ stochastic simulation algorithms and Brownian dynamics, respectively. Through theoretical analysis, we have developed an algorithm to identify if the system is reaction-controlled, diffusion-controlled or is in an intermediate regime. The time-step size is chosen accordingly at each step of the simulation. We have used three examples to demonstrate the accuracy and robustness of the proposed algorithm. The first example deals with diffusion of two chemical species undergoing an irreversible bimolecular reaction. It is used to validate our algorithm by comparing its results with the solution obtained from a corresponding deterministic partial differential equation at low and high number of molecules. In this example, we also compare the results from our method to those obtained using a Gillespie multi-particle (GMP) method. The second example, which models simplified RNA synthesis, is used to study the performance of our algorithm in reaction- and diffusion-controlled regimes and to investigate the effects of local inhomogeneity. The third example models reaction-diffusion of CheY molecules through the cytoplasm of Escherichia coli during chemotaxis. It is used to compare the algorithm's performance against the GMP method. Our analysis demonstrates that the proposed algorithm enables accurate simulation of the kinetics of complex and spatially heterogeneous systems. It is also computationally more efficient than commonly used alternatives, such as the GMP method.

Autonomous learning of nonlocal stochastic neuron dynamics.

Neuronal dynamics is driven by externally imposed or internally generated random excitations/noise, and is often described by systems of random or stochastic ordinary differential equations. Such systems admit a distribution of solutions, which is (partially) characterized by the single-time joint probability density function (PDF) of system states. It can be used to calculate such information-theoretic quantities as the mutual information between the stochastic stimulus and various internal states of the neuron (e.g., membrane potential), as well as various spiking statistics. When random excitations are modeled as Gaussian white noise, the joint PDF of neuron states satisfies exactly a Fokker-Planck equation. However, most biologically plausible noise sources are correlated (colored). In this case, the resulting PDF equations require a closure approximation. We propose two methods for closing such equations: a modified nonlocal large-eddy-diffusivity closure and a data-driven closure relying on sparse regression to learn relevant features. The closures are tested for the stochastic non-spiking leaky integrate-and-fire and FitzHugh-Nagumo (FHN) neurons driven by sine-Wiener noise. Mutual information and total correlation between the random stimulus and the internal states of the neuron are calculated for the FHN neuron.