Search

We develop computational methods for the simulation of osmotic swelling phenomena relevant to microscopic vesicles containing transformable solute molecules. We introduce stochastic immersed boundary methods (SIBMs) that can capture osmotically driven fluid transport through semipermeable elastic membranes subject to thermal fluctuations. We also develop numerical methods to handle within SIBMs an elastic shell model for a neo-Hookean material. Our extended SIBM allows for capturing osmotic swelling phenomena driven by concentration changes and interactions between a discrete collection of confined particles while accounting for the thermal fluctuations of the semipermeable membrane and the hydrodynamic transport of solvent. We use our computational methods to investigate osmotic phenomena in regimes that go beyond the classical Van't Hoff theory. We develop statistical mechanics theories for osmotic swelling of vesicles when there are significant interactions between particles that can transform over time. We validate our theoretical results against detailed computational simulations. Our methods are expected to be useful for a wide class of applications allowing for the simulation of osmotically driven flows, thermally fluctuating semipermeable elastic structures, and solute interactions.

In this work, we examine how volume exclusion caused by regions of high chromatin density might influence the time required for proteins to find specific DNA binding sites. The spatial variation of chromatin density within mouse olfactory sensory neurons i

Cardiac fluid dynamics fundamentally involves interactions between complex blood flows and the structural deformations of the muscular heart walls and the thin valve leaflets. There has been longstanding scientific, engineering, and medical interest in creating mathematical models of the heart that capture, explain, and predict these fluid-structure interactions (FSIs). However, existing computational models that account for interactions among the blood, the actively contracting myocardium, and the valves are limited in their abilities to predict valve performance, capture fine-scale flow features, or use realistic descriptions of tissue biomechanics. Here we introduce and benchmark a comprehensive mathematical model of cardiac FSI in the human heart. A unique feature of our model is that it incorporates biomechanically detailed descriptions of all major cardiac structures that are calibrated using tensile tests of human tissue specimens to reflect the heart's microstructure. Further, it is the first FSI model of the heart that provides anatomically and physiologically detailed representations of all four cardiac valves. We demonstrate that this integrative model generates physiologic dynamics, including realistic pressure-volume loops that automatically capture isovolumetric contraction and relaxation, and that its responses to changes in loading conditions are consistent with the Frank-Starling mechanism. These complex relationships emerge intrinsically from interactions within our comprehensive description of cardiac physiology. Such models can serve as tools for predicting the impacts of medical interventions. They also can provide platforms for mechanistic studies of cardiac pathophysiology and dysfunction, including congenital defects, cardiomyopathies, and heart failure, that are difficult or impossible to perform in patients., (© The Author(s) 2024. Published by Oxford University Press on behalf of National Academy of Sciences.)

Synaptic plasticity [long-term potentiation/depression (LTP/D)], is a cellular mechanism underlying learning. Two distinct types of early LTP/D (E-LTP/D), acting on very different time scales, have been observed experimentally-spike timing dependent plasticity (STDP), on time scales of tens of ms; and behavioral time scale synaptic plasticity (BTSP), on time scales of seconds. BTSP is a candidate for a mechanism underlying rapid learning of spatial location by place cells. Here, a computational model of the induction of E-LTP/D at a spine head of a synapse of a hippocampal pyramidal neuron is developed. The single-compartment model represents two interacting biochemical pathways for the activation (phosphorylation) of the kinase (CaMKII) with a phosphatase, with ion inflow through channels (NMDAR, CaV1,Na). The biochemical reactions are represented by a deterministic system of differential equations, with a detailed description of the activation of CaMKII that includes the opening of the compact state of CaMKII. This single model captures realistic responses (temporal profiles with the differing timescales) of STDP and BTSP and their asymmetries. The simulations distinguish several mechanisms underlying STDP vs. BTSP, including i) the flow of [Formula: see text] through NMDAR vs. CaV1 channels, and ii) the origin of several time scales in the activation of CaMKII. The model also realizes a priming mechanism for E-LTP that is induced by [Formula: see text] flow through CaV1.3 channels. Once in the spine head, this small additional [Formula: see text] opens the compact state of CaMKII, placing CaMKII ready for subsequent induction of LTP., Competing Interests: Competing interests statement:The authors declare no competing interest.

This paper describes computer models of three interventions used for treating refractory pulmonary hypertension (RPH). These procedures create either an atrial septal defect, a ventricular septal defect or, in the case of a Potts shunt, a patent ductus arteriosus. The aim in all three cases is to generate a right-to-left shunt, allowing for either pressure or volume unloading of the right side of the heart in the setting of right ventricular failure, while maintaining cardiac output. These shunts are created, however, at the expense of introducing de-oxygenated blood into the systemic circulation, thereby lowering the systemic arterial oxygen saturation. The models developed in this paper are based on compartmental descriptions of human hemodynamics and oxygen transport. An important parameter included in our models is the cross-sectional area of the surgically created defect. Numerical simulations are performed to compare different interventions and various shunt sizes and to assess their impact on hemodynamic variables and oxygen saturations. We also create a model for exercise and use it to study exercise tolerance in simulated pre-intervention and post-intervention RPH patients., (© The Author(s) 2022. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.)

In this paper, we develop a pulsatile compartmental model of the Fontan circulation and use it to explore the effects of a fenestration added to this physiology. A fenestration is a shunt between the systemic and pulmonary veins that is added either at the time of Fontan conversion or at a later time for the treatment of complications. This shunt increases cardiac output and decreases systemic venous pressure. However, these hemodynamic benefits are achieved at the expense of a decrease in the arterial oxygen saturation. The model developed in this paper incorporates fenestration size as a parameter and describes both blood flow and oxygen transport. It is calibrated to clinical data from Fontan patients, and we use it to study the impact of a fenestration on several hemodynamic variables, including systemic oxygen availability, effective oxygen availability, and systemic venous pressure. In certain scenarios corresponding to high-risk Fontan physiology, we demonstrate the existence of a range of fenestration sizes in which the systemic oxygen availability remains relatively constant while the systemic venous pressure decreases., Competing Interests: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest., (Copyright © 2022 Ahmad, Jin, Penny, Rusin, Peskin and Puelz.)

Based on von Neumann's model of an economy characterized by processes and goods, we add to that model a component representing capital equipment. We assume that the need for capital equipment by any process is proportional to the rate at which that process is running, and therefore an increase in rate requires that capital equipment be purchased, whereas a decrease in rate allows capital equipment to be sold. We thereby construct a continuous-time dynamical model, which we use to investigate the evolution of economic diversity under two price equilibrium scenarios: the first with non-negative prices and non-positive excess demands; the second with enforced market clearing and with prices allowed to be negative. The second scenario represents an economy in which recycling is required, so that excess supply cannot be discarded. We prove that at any time during the progression of the model economy, the solution to each of the two price equilibrium problems exists, and that non-uniqueness of the solution, if any, does not affect the development of the model economy. We compare matched model economies under the two scenarios by simulating their respective evolutions. In each case, the model economy experiences a process of selection and matures to a state of balanced growth, with a higher growth rate when excess supply is discarded, but with greater economic diversity with enforced recycling. The robustness of these qualitative results is demonstrated by repeated trials of simulations on matched pairs of model economies with different randomly chosen parameters., Competing Interests: The authors have declared that no competing interests exist.

This paper focuses on the derivation and simulation of mathematical models describing new plasma fraction in blood for patients undergoing simultaneous extracorporeal membrane oxygenation and therapeutic plasma exchange. Models for plasma exchange with either veno-arterial or veno-venous extracorporeal membrane oxygenation are considered. Two classes of models are derived for each case, one in the form of an algebraic delay equation and another in the form of a system of delay differential equations. In special cases, our models reduce to single compartment ones for plasma exchange that have been validated with experimental data (Randerson et al., 1982, Artif. Organs, 6, 43-49). We also show that the algebraic differential equations are forward Euler discretizations of the delay differential equations, with timesteps equal to transit times through model compartments. Numerical simulations are performed to compare different model types, to investigate the impact of plasma device port switching on the efficiency of the exchange process, and to study the sensitivity of the models to their parameters., (© The Author(s) 2021. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.)

Vascular access connection configurations during tandem extracorporeal membrane oxygenation (ECMO) and therapeutic plasma exchange (TPE) may impact exchange kinetics. In these tandem procedures, typically the TPE inlet line is proximal to the TPE return line with respect to blood flow in the ECMO device, maximizing the opportunity for replacement fluid homogenization within the ECMO circuit. However, if TPE inlet and return line connections are switched, recirculation-a phenomenon in which replacement fluid leaving the TPE return line is prematurely drawn into the TPE inlet line prior to satisfactory homogenization within the ECMO circuit-will occur. Such recirculation could diminish TPE efficacy in patients on ECMO and mitigate therapeutic benefits. Using a mathematical model of recirculation in tandem ECMO and TPE, we demonstrate that the predicted impact of recirculation is negligible and vascular access connection positioning does not appear to be a point of clinical concern with regard to TPE kinetics., (© 2020 Wiley Periodicals LLC.)

For a chemical signal to propagate across a cell, it must navigate a tortuous environment involving a variety of organelle barriers. In this work we study mathematical models for a basic chemical signal, the arrival times at the nuclear membrane of proteins that are activated at the cell membrane and diffuse throughout the cytosol. Organelle surfaces within human B cells are reconstructed from soft X-ray tomographic images, and modeled as reflecting barriers to the molecules' diffusion. We show that signal inactivation sharpens signals, reducing variability in the arrival time at the nuclear membrane. Inactivation can also compensate for an observed slowdown in signal propagation induced by the presence of organelle barriers, leading to arrival times at the nuclear membrane that are comparable to models in which the cytosol is treated as an open, empty region. In the limit of strong signal inactivation this is achieved by filtering out molecules that traverse non-geodesic paths., Competing Interests: The authors have declared that no competing interests exist.

Widespread use of antibiotics has resulted in an increase in antimicrobial-resistant microorganisms. Although not all bacterial contact results in infection, patients can become asymptomatically colonized, increasing the risk of infection and pathogen transmission. Consequently, many institutions have begun active surveillance, but in non-research settings, the resulting data are often incomplete and may include non-random testing, making conventional epidemiological analysis problematic. We describe a mathematical model and inference method for in-hospital bacterial colonization and transmission of carbapenem-resistant Enterobacteriaceae that is tailored for analysis of active surveillance data with incomplete observations. The model and inference method make use of the full detailed state of the hospital unit, which takes into account the colonization status of each individual in the unit and not only the number of colonized patients at any given time. The inference method computes the exact likelihood of all possible histories consistent with partial observations (despite the exponential increase in possible states that can make likelihood calculation intractable for large hospital units), includes techniques to improve computational efficiency, is tested by computer simulation, and is applied to active surveillance data from a 13-bed rehabilitation unit in New York City. The inference method for exact likelihood calculation is applicable to other Markov models incorporating incomplete observations. The parameters that we identify are the patient-patient transmission rate, pre-existing colonization probability, and prior-to-new-patient transmission probability. Besides identifying the parameters, we predict the effects on the total prevalence (0.07 of the total colonized patient-days) of changing the parameters and estimate the increase in total prevalence attributable to patient-patient transmission (0.02) above the baseline pre-existing colonization (0.05). Simulations with a colonized versus uncolonized long-stay patient had 44% higher total prevalence, suggesting that the long-stay patient may have been a reservoir of transmission. High-priority interventions may include isolation of incoming colonized patients and repeated screening of long-stay patients., Competing Interests: The authors have declared that no competing interests exist.

This work is concerned with modeling and simulation of the mitral valve, one of the four valves in the human heart. The valve is composed of leaflets, the free edges of which are supported by a system of chordae, which themselves are anchored to the papillary muscles inside the left ventricle. First, we examine valve anatomy and present the results of original dissections. These display the gross anatomy and information on fiber structure of the mitral valve. Next, we build a model valve following a design-based methodology, meaning that we derive the model geometry and the forces that are needed to support a given load and construct the model accordingly. We incorporate information from the dissections to specify the fiber topology of this model. We assume the valve achieves mechanical equilibrium while supporting a static pressure load. The solution to the resulting differential equations determines the pressurized configuration of the valve model. To complete the model, we then specify a constitutive law based on a stress-strain relation consistent with experimental data that achieves the necessary forces computed in previous steps. Finally, using the immersed boundary method, we simulate the model valve in fluid in a computer test chamber. The model opens easily and closes without leak when driven by physiological pressures over multiple beats. Further, its closure is robust to driving pressures that lack atrial systole or are much lower or higher than normal., (© 2019 John Wiley & Sons, Ltd.)

In this paper, we consider a stochastic molecular circadian oscillator described by a sequence of biological reactions and its deterministic kinetics governed by a system of ordinary differential equations in the limit of large numbers of molecules. The oscillations in the model are generated by negative feedback regulation of a gene. The focus of this paper is the entrainment of the oscillator by a periodic light signal that affects the maximal transcription rate of the gene. We introduce two scalings of the model parameters that provide independent control over the natural frequency of the oscillator and the relative noise level. We study entrainment in two ways: by visualizing the stochastic limit cycle in various projections of the discrete phase space of the system and by evaluating the maximum of the normalized cross correlation of the light signal with the number of protein molecules in the cell. The visualization method ignores the phase of the oscillator, and we find in this way that entrainment has a subtle organizing effect on the limit cycle as a whole. The cross correlation results reveal an interval of natural frequencies of the oscillator surrounding the frequency of the light signal within which maximal entrainment occurs with rather sharp drops in entrainment at the edges of this interval. The width of the interval of maximal entrainment increases with the amplitude of the light signal. These statements are applicable both to the stochastic oscillator and to its deterministic limit, but the results are most clear-cut in the deterministic case and degrade from there as the relative noise level increases.

The exact mechanism to orchestrate the action of hundreds of dynein motor proteins to generate wave-like ciliary beating remains puzzling and has fascinated many scientists. We present a 3D model of a cilium and the simulation of its beating in a fluid environment. The model cilium obeys a simple geometric constraint that arises naturally from the microscopic structure of a real cilium. This constraint allows us to determine the whole 3D structure at any instant in terms of the configuration of a single space curve. The tensions of active links, which model the dynein motor proteins, follow a postulated dynamical law, and together with the passive elasticity of microtubules, this dynamical law is responsible for the ciliary motions. In particular, our postulated tension dynamics lead to the instability of a symmetrical steady state, in which the cilium is straight and its active links are under equal tensions. The result of this instability is a stable, wave-like, limit cycle oscillation. We have also investigated the fluid-structure interaction of cilia using the immersed boundary (IB) method. In this setting, we see not only coordination within a single cilium but also, coordinated motion, in which multiple cilia in an array organize their beating to pump fluid, in particular by breaking phase synchronization., Competing Interests: The authors declare no conflict of interest.

The Immersed Boundary (IB) method is a mathematical framework for constructing robust numerical methods to study fluid-structure interaction in problems involving an elastic structure immersed in a viscous fluid. The IB formulation uses an Eulerian representation of the fluid and a Lagrangian representation of the structure. The Lagrangian and Eulerian frames are coupled by integral transforms with delta function kernels. The discretized IB equations use approximations to these transforms with regularized delta function kernels to interpolate the fluid velocity to the structure, and to spread structural forces to the fluid. It is well-known that the conventional IB method can suffer from poor volume conservation since the interpolated Lagrangian velocity field is not generally divergence-free, and so this can cause spurious volume changes. In practice, the lack of volume conservation is especially pronounced for cases where there are large pressure differences across thin structural boundaries. The aim of this paper is to greatly reduce the volume error of the IB method by introducing velocity-interpolation and force-spreading schemes with the properties that the interpolated velocity field in which the structure moves is at least C 1 and satisfies a continuous divergence-free condition, and that the force-spreading operator is the adjoint of the velocity-interpolation operator. We confirm through numerical experiments in two and three spatial dimensions that this new IB method is able to achieve substantial improvement in volume conservation compared to other existing IB methods, at the expense of a modest increase in the computational cost. Further, the new method provides smoother Lagrangian forces (tractions) than traditional IB methods. The method presented here is restricted to periodic computational domains. Its generalization to non-periodic domains is important future work.

We simulate deformable red blood cells in the microcirculation using the immersed boundary method with a cytoskeletal model that incorporates structural details revealed by tomographic images. The elasticity of red blood cells is known to be supplied by both their lipid bilayer membranes, which resist bending and local changes in area, and their cytoskeletons, which resist in-plane shear. The cytoskeleton consists of spectrin tetramers that are tethered to the lipid bilayer by ankyrin and by actin-based junctional complexes. We model the cytoskeleton as a random geometric graph, with nodes corresponding to junctional complexes and with edges corresponding to spectrin tetramers such that the edge lengths are given by the end-to-end distances between nodes. The statistical properties of this graph are based on distributions gathered from three-dimensional tomographic images of the cytoskeleton by a segmentation algorithm. We show that the elastic response of our model cytoskeleton, in which the spectrin polymers are treated as entropic springs, is in good agreement with the experimentally measured shear modulus. By simulating red blood cells in flow with the immersed boundary method, we compare this discrete cytoskeletal model to an existing continuum model and predict the extent to which dynamic spectrin network connectivity can protect against failure in the case of a red cell subjected to an applied strain. The methods presented here could form the basis of disease- and patient-specific computational studies of hereditary diseases affecting the red cell cytoskeleton.

Neuroendocrine control of reproduction by brain-secreted pulses of gonadotropin-releasing hormone (GnRH) represents a longstanding puzzle about extracellular signal decoding mechanisms. GnRH regulates the pituitary gonadotropin's follicle-stimulating hormone (FSH) and luteinizing hormone (LH), both of which are heterodimers specified by unique β subunits (FSHβ/LHβ). Contrary to Lhb , Fshb gene induction has a preference for low-frequency GnRH pulses. To clarify the underlying regulatory mechanisms, we developed three biologically anchored mathematical models: 1) parallel activation of Fshb inhibitory factors ( e.g. inhibin α and VGF nerve growth factor-inducible), 2) activation of a signaling component with a refractory period ( e.g. G protein), and 3) inactivation of a factor needed for Fshb induction ( e.g. growth differentiation factor 9). Simulations with all three models recapitulated the Fshb expression levels obtained in pituitary gonadotrope cells perifused with varying GnRH pulse frequencies. Notably, simulations altering average concentration, pulse duration, and pulse frequency revealed that the apparent frequency-dependent pattern of Fshb expression in model 1 actually resulted from variations in average GnRH concentration. In contrast, models 2 and 3 showed "true" pulse frequency sensing. To resolve which components of this GnRH signal induce Fshb , we developed a high-throughput parallel experimental system. We analyzed over 4,000 samples in experiments with varying near-physiological GnRH concentrations and pulse patterns. Whereas Egr1 and Fos genes responded only to variations in average GnRH concentration, Fshb levels were sensitive to both average concentration and true pulse frequency. These results provide a foundation for understanding the role of multiple regulatory factors in modulating Fshb gene activity., (© 2017 by The American Society for Biochemistry and Molecular Biology, Inc.)

The helical flagella that are attached to the cell body of bacteria such as Escherichia coli and Salmonella typhimurium allow the cell to swim in a fluid environment. These flagella are capable of polymorphic transformation in that they take on various helical shapes that differ in helical pitch, radius, and chirality. We present a mathematical model of a single flagellum described by Kirchhoff rod theory that is immersed in a fluid governed by Stokes equations. We perform numerical simulations to demonstrate two mechanisms by which polymorphic transformation can occur, as observed in experiments. First, we consider a flagellar filament attached to a rotary motor in which transformations are triggered by a reversal of the direction of motor rotation [L. Turner et al., J. Bacteriol. 182, 2793 (2000)10.1128/JB.182.10.2793-2801.2000]. We then consider a filament that is fixed on one end and immersed in an external fluid flow [H. Hotani, J. Mol. Biol. 156, 791 (1982)10.1016/0022-2836(82)90142-5]. The detailed dynamics of the helical flagellum interacting with a viscous fluid is discussed and comparisons with experimental and theoretical results are provided.

Background: Echocardiography-derived linear dimensions offer straightforward indices of right ventricular (RV) structure but have not been systematically compared with RV volumes on cardiac magnetic resonance (CMR)., Methods: Echocardiography and CMR were interpreted among patients with coronary artery disease imaged via prospective (90%) and retrospective (10%) registries. For echocardiography, American Society of Echocardiography-recommended RV dimensions were measured in apical four-chamber (basal RV width, mid RV width, and RV length), parasternal long-axis (proximal RV outflow tract [RVOT]), and short-axis (distal RVOT) views. For CMR, RV end-diastolic volume and RV end-systolic volume were quantified using border planimetry., Results: Two hundred seventy-two patients underwent echocardiography and CMR within a narrow interval (0.4 ± 1.0 days); complete acquisition of all American Society of Echocardiography-recommended dimensions was feasible in 98%. All echocardiographic dimensions differed between patients with and those without RV dilation on CMR (P < .05). Basal RV width (r = 0.70), proximal RVOT width (r = 0.68), and RV length (r = 0.61) yielded the highest correlations with RV end-diastolic volume on CMR; end-systolic dimensions yielded similar correlations (r = 0.68, r = 0.66, and r = 0.65, respectively). In multivariate regression, basal RV width (regression coefficient = 1.96 per mm; 95% CI, 1.22-2.70; P < .001), RV length (regression coefficient = 0.97; 95% CI, 0.56-1.37; P < .001), and proximal RVOT width (regression coefficient = 2.62; 95% CI, 1.79-3.44; P < .001) were independently associated with CMR RV end-diastolic volume (r = 0.80). RV end-systolic volume was similarly associated with echocardiographic dimensions (basal RV width: 1.59 per mm [95% CI, 1.06-2.13], P < .001; RV length: 1.00 [95% CI, 0.66-1.34], P < .001; proximal RVOT width: 1.80 [95% CI, 1.22-2.39], P < .001) (r = 0.79)., Conclusions: RV linear dimensions provide readily obtainable markers of RV chamber size. Proximal RVOT and basal width are independently associated with CMR volumes, supporting the use of multiple linear dimensions when assessing RV size on echocardiography., Competing Interests: Conflicts of Interests Disclosure: None, (Copyright © 2016 American Society of Echocardiography. Published by Elsevier Inc. All rights reserved.)

Determining the cellular basis of brain growth is an important problem in developmental neurobiology. In the mammalian brain, the cerebellum is particularly amenable to studies of growth because it contains only a few cell types, including the granule cells, which are the most numerous neuronal subtype. Furthermore, in the mouse cerebellum granule cells are generated from granule cell precursors (gcps) in the external granule layer (EGL), from 1 day before birth until about 2 weeks of age. The complexity of the underlying cellular processes (multiple cell behaviors, three spatial dimensions, time-dependent changes) requires a quantitative framework to be fully understood. In this paper, a differential equation-based model is presented, which can be used to estimate temporal changes in granule cell numbers in the EGL. The model includes the proliferation of gcps and their differentiation into granule cells, as well as the process by which granule cells leave the EGL. Parameters describing these biological processes were derived from fitting the model to histological data. This mathematical model should be useful for understanding altered gcp and granule cell behaviors in mouse mutants with abnormal cerebellar development and cerebellar cancers.

The probabilistic nature of neurotransmitter release in synapses is believed to be one of the most significant sources of noise in the central nervous system. We show how p0, the probability of release per docked vesicle when an action potential arrives, affects the dynamics of the rate of vesicle release in response to changes in the rate of arrival of action potentials. Furthermore, we examine the theoretical capability of a synapse in the estimation of desired signals using information from the stochastic vesicle release events under the framework of optimal linear filter theory. We find that a small p0, such as 0.1, reduces the error in the reconstruction of the input, or in the reconstruction of the time derivative of the input, from the time series of vesicle release events. Our results imply that the probabilistic nature of synaptic vesicle release plays a direct functional role in synaptic transmission.

In this work, we consider a stochastic oscillator described by a discrete-state continuous-time Markov chain, in which the states are arranged in a circle, and there is a constant probability per unit time of jumping from one state to the next in a specified direction around the circle. At each of a sequence of equally spaced times, the oscillator has a specified probability of being reset to a particular state. The focus of this work is the entrainment of the oscillator by this periodic but stochastic stimulus. We consider a distinguished limit, in which (i) the number of states of the oscillator approaches infinity, as does the probability per unit time of jumping from one state to the next, so that the natural mean period of the oscillator remains constant, (ii) the resetting probability approaches zero, and (iii) the period of the resetting signal approaches a multiple, by a ratio of small integers, of the natural mean period of the oscillator. In this distinguished limit, we use analytic and numerical methods to study the extent to which entrainment occurs.

A novel multiscale mathematical and computational model of the pulmonary circulation is presented and used to analyse both arterial and venous pressure and flow. This work is a major advance over previous studies by Olufsen et al. (Ann Biomed Eng 28:1281-1299, 2012) which only considered the arterial circulation. For the first three generations of vessels within the pulmonary circulation, geometry is specified from patient-specific measurements obtained using magnetic resonance imaging (MRI). Blood flow and pressure in the larger arteries and veins are predicted using a nonlinear, cross-sectional-area-averaged system of equations for a Newtonian fluid in an elastic tube. Inflow into the main pulmonary artery is obtained from MRI measurements, while pressure entering the left atrium from the main pulmonary vein is kept constant at the normal mean value of 2 mmHg. Each terminal vessel in the network of 'large' arteries is connected to its corresponding terminal vein via a network of vessels representing the vascular bed of smaller arteries and veins. We develop and implement an algorithm to calculate the admittance of each vascular bed, using bifurcating structured trees and recursion. The structured-tree models take into account the geometry and material properties of the 'smaller' arteries and veins of radii ≥ 50 μm. We study the effects on flow and pressure associated with three classes of pulmonary hypertension expressed via stiffening of larger and smaller vessels, and vascular rarefaction. The results of simulating these pathological conditions are in agreement with clinical observations, showing that the model has potential for assisting with diagnosis and treatment for circulatory diseases within the lung.

Membrane current through voltage-sensitive calcium ion channels at the postsynaptic density of a dendritic spine is investigated. To simulate the ion channels that carry such current and the resulting temporal and spatial distribution of concentration, current, and voltage within the dendritic spine, the immersed boundary method with electrodiffusion is applied. In this simulation method a spatially continuous chemical potential barrier is used to simulate the influence of the membrane on each species of ion. The amplitudes of these barriers can be regulated to simulate channel gating. Here we introduce this methodology in a one-dimensional setting. First, we study the current-voltage relationship obtained with fixed chemical potential barriers. Next, we simulate stochastic ion-channel gating in a calcium channel with multiple subunits, and observe the diffusive wave of calcium entry within the dendritic spine that follows channel opening. This work lays the foundation for future three-dimensional studies of electrodiffusion and advection electrodiffusion in dendritic spines., (© 2013 Elsevier Ltd. All rights reserved.)

Electrical signaling is a fast mode of communication for cells within an organism. We are concerned here with the formulation and analysis of mathematical models that are used to describe this important class of physiological processes. These models generally take the form of partial differential equations that are descendants of those introduced by Hodgkin and Huxley to describe the propagation of an action potential along the squid giant axon. We review that work here and then go on to describe more recent variations on the Hodgkin-Huxley theme, including the three-dimensional bidomain (and monodomain) equations for cardiac electrophysiology, multiscale models for the heart that take cellular structure into account near the action potential wavefront, and finally a more detailed reformulation of electrophysiology in terms of electrodiffusion.

In this work, we examine how volume exclusion caused by regions of high chromatin density might influence the time required for proteins to find specific DNA binding sites. The spatial variation of chromatin density within mouse olfactory sensory neurons is determined from soft X-ray tomography reconstructions of five nuclei. We show that there is a division of the nuclear space into regions of low-density euchromatin and high-density heterochromatin. Volume exclusion experienced by a diffusing protein caused by this varying density of chromatin is modeled by a repulsive potential. The value of the potential at a given point in space is chosen to be proportional to the density of chromatin at that location. The constant of proportionality, called the volume exclusivity, provides a model parameter that determines the strength of volume exclusion. Numerical simulations demonstrate that the mean time for a protein to locate a binding site localized in euchromatin is minimized for a finite, nonzero volume exclusivity. For binding sites in heterochromatin, the mean time is minimized when the volume exclusivity is zero (the protein experiences no volume exclusion). An analytical theory is developed to explain these results. The theory suggests that for binding sites in euchromatin there is an optimal level of volume exclusivity that balances a reduction in the volume searched in finding the binding site, with the height of effective potential barriers the protein must cross during the search process.

We present a generalization of a population density approach for modeling and analysis of stochastic gene expression. In the model, the gene of interest fluctuates stochastically between an inactive state, in which transcription cannot occur, and an active state, in which discrete transcription events occur; and the individual mRNA molecules are degraded stochastically in an independent manner. This sort of model in simplest form with exponential dwell times has been used to explain experimental estimates of the discrete distribution of random mRNA copy number. In our generalization, the random dwell times in the inactive and active states, T&#95;{0} and T&#95;{1}, respectively, are independent random variables drawn from any specified distributions. Consequently, the probability per unit time of switching out of a state depends on the time since entering that state. Our method exploits a connection between the fully discrete random process and a related continuous process. We present numerical methods for computing steady-state mRNA distributions and an analytical derivation of the mRNA autocovariance function. We find that empirical estimates of the steady-state mRNA probability mass function from Monte Carlo simulations of laboratory data do not allow one to distinguish between underlying models with exponential and nonexponential dwell times in some relevant parameter regimes. However, in these parameter regimes and where the autocovariance function has negative lobes, the autocovariance function disambiguates the two types of models. Our results strongly suggest that temporal data beyond the autocovariance function is required in general to characterize gene switching.

We describe and analyze a model for a stochastic pulse-coupled neuronal network with many sources of randomness: random external input, potential synaptic failure, and random connectivity topologies. We show that different classes of network topologies give rise to qualitatively different types of synchrony: uniform (Erdős-Rényi) and "small-world" networks give rise to synchronization phenomena similar to that in "all-to-all" networks (in which there is a sharp onset of synchrony as coupling is increased); in contrast, in "scale-free" networks the dependence of synchrony on coupling strength is smoother. Moreover, we show that in the uniform and small-world cases, the fine details of the network are not important in determining the synchronization properties; this depends only on the mean connectivity. In contrast, for scale-free networks, the dynamics are significantly affected by the fine details of the network; in particular, they are significantly affected by the local neighborhoods of the "hubs" in the network.

Flagellar bundling is an important aspect of locomotion in bacteria such as Escherichia coli. To study the hydrodynamic behavior of helical flagella, we present a computational model that is based on the geometry of the bacterial flagellar filament at the micrometer scale. We consider two model flagella, each of which has a rotary motor at its base with the rotation rate of the motor set at 100 Hz. Bundling occurs when both flagella are left-handed helices turning counterclockwise (when viewed from the nonmotor end of the flagellum looking back toward the motor) or when both flagella are right-handed helices turning clockwise. Helical flagella of the other combinations of handedness and rotation direction do not bundle.

Within the nuclei of eukaryotic cells, the density of chromatin is nonuniform. We study the influence of this nonuniform density, which we derive from microscopic images [Schermelleh L, et al. (2008) Science 320:1332-1336], on the diffusion of proteins within the nucleus, under the hypothesis that chromatin density is proportional to an effective potential that tends to exclude the diffusing protein from regions of high chromatin density. The constant of proportionality, which we call the volume exclusivity of chromatin, is a model parameter that we can tune to study the influence of such volume exclusivity on the random time required for a diffusing particle to find its target. We consider randomly chosen binding sites located in regions of low (20th-30th percentile) chromatin density, and we compute the median time to find such a binding site by a protein that enters the nucleus at a randomly chosen nuclear pore. As the volume exclusivity of chromatin increases from zero, we find that the median time needed to reach the target binding site at first decreases to a minimum, and then increases again as the volume exclusivity of chromatin increases further. Random permutation of the voxel values of chromatin density abolishes the minimum, thus demonstrating that the speedup seen with increasing volume exclusivity at low to moderate volume exclusivity is dependent upon the spatial structure of chromatin within the nucleus.

We derive a homogenized description of the electrical communication along a single strand of myocytes in the presence of gap-junctional and electric-field coupling. In the model, cells are electrically coupled through narrow clefts that are resistively connected to extracellular space. Cells are also coupled directly through gap junctions. We perform numerical simulations of this full model and its homogenization. We observe that the full and homogenized descriptions agree when gap-junctional coupling is at physiologically normal levels. When gap-junctional coupling is low, the two descriptions disagree. In this case, only the full model captures the behavior that the ephaptic mechanism can speed up action potential propagation. A strength of our homogenized description is that it is a macroscale model that can account for the preferential localization of Na+ channels at the ends of cells.

Objective: Cerebral autoregulation maintains a relatively stable cerebral blood flow over a range of perfusion pressures. During exercise, regional cerebral blood flow may be elevated in particular areas of the brain. This case report presents the impact of aerobic exercise on intracranially measured pressure and brain tissue oxygenation in an adult human., Clincial Presentation: A 30-year-old man with idiopathic intracranial hypertension treated with cerebrospinal fluid diversion was monitored with a Licox intracranial brain oxygen and pressure monitor (Integra NeuroSciences Corporation, Plainsboro, New Jersey) for refractory nonpostural headaches exacerbated after exercise. He performed trials of running and bicycling to provoke his headaches. The patient's mean intracranial pressure remained stable during exercise despite elevated cerebral perfusion pressures. Regional cerebral oxygen tension levels were strictly regulated to a level of approximately 39 mm Hg during steady state aerobic exercise, with transient increases up to 90 mm Hg at the onset and termination of activity., Conclusion: Our results suggest that cerebral autoregulation appears to maintain constant cerebral oxygen tension during exercise. Further, we note transient cerebral hyperoxygenation at the onset of exercise as autoregulation "turns on" and at the termination of exercise. We present a quantitative interpretation of the post-exercise hyperoxygenation phase based on Fick's principle. We are the first to demonstrate cortical hyperoxygenation in a human breathing natural air without oxygen supplementation.

We describe an immersed boundary method for problems of fluid-solute-structure interaction. The numerical scheme employs linearly implicit timestepping, allowing for the stable use of timesteps that are substantially larger than those permitted by an explicit method, and local mesh refinement, making it feasible to resolve the steep gradients associated with the space charge layers as well as the chemical potential, which is used in our formulation to control the permeability of the membrane to the (possibly charged) solute. Low Reynolds number fluid dynamics are described by the time-dependent incompressible Stokes equations, which are solved by a cell-centered approximate projection method. The dynamics of the chemical species are governed by the advection-electrodiffusion equations, and our semi-implicit treatment of these equations results in a linear system which we solve by GMRES preconditioned via a fast adaptive composite-grid (FAC) solver. Numerical examples demonstrate the capabilities of this methodology, as well as its convergence properties.

In normal somatic cells, telomere length shortens with each cell replication. This progressive shortening is associated with cellular senescence and apoptosis. Germ cells, stem cells, and the majority of cancer cells express telomerase, an enzyme that extends telomere length and, when expressed at sufficient levels, can immortalize or extend the life span of a cell line. It is believed that telomeres switch between two states: capped and uncapped. The telomere state determines its accessibility to telomerase and also the onset of senescence. One hypothesis is that the t loop, a large lariat-like structure, represents the capped state. In this paper we model a telomere state on the basis of the biophysics of t-loop formation, allowing us to develop a single mathematical model that accounts for two processes: telomere length regulation for telomerase positive cells and cellular senescence in somatic cells. The model predicts the steady-state length distribution for telomerase positive cells, describes the time evolution of telomere length, and computes the life span of a cell line on the basis of the levels of TRF2 and telomerase expression. The model reproduces a wide range of experimental behavior and fits data from immortal cell lines (HeLa S3 and 293T) and somatic cells (human diploid fibroblasts) well. We conclude that the t loop as the capped state is a quantitatively reasonable model of telomere length regulation and cellular senescence.

In this paper we study the error rate of RNA synthesis in the look-ahead model for the random walk of RNA polymerase along DNA during transcription. The model's central assumption is the existence of a window of activity in which ribonucleoside triphosphates (rNTPs) bind reversibly to the template DNA strand before being hydrolyzed and linked covalently to the nascent RNA chain. An unknown, but important, integer parameter of this model is the window size w. Here, we use mathematical analysis and computer simulation to study the rate at which transcriptional errors occur as a function of w. We find dramatic reduction in the error rate of transcription as w increases, especially for small values of w. The error reduction method provided by look-ahead occurs before hydrolysis and covalent linkage of rNTP to the nascent RNA chain, and is therefore distinct from error correction mechanisms that have previously been considered.

Many physiological responses elicited by neuronal spikes-intracellular calcium transients, synaptic potentials, muscle contractions-are built up of discrete, elementary responses to each spike. However, the spikes occur in trains of arbitrary temporal complexity, and each elementary response not only sums with previous ones, but can itself be modified by the previous history of the activity. A basic goal in system identification is to characterize the spike-response transform in terms of a small number of functions-the elementary response kernel and additional kernels or functions that describe the dependence on previous history-that will predict the response to any arbitrary spike train. Here we do this by developing further and generalizing the "synaptic decoding" approach of Sen et al. (1996). Given the spike times in a train and the observed overall response, we use least-squares minimization to construct the best estimated response and at the same time best estimates of the elementary response kernel and the other functions that characterize the spike-response transform. We avoid the need for any specific initial assumptions about these functions by using techniques of mathematical analysis and linear algebra that allow us to solve simultaneously for all of the numerical function values treated as independent parameters. The functions are such that they may be interpreted mechanistically. We examine the performance of the method as applied to synthetic data. We then use the method to decode real synaptic and muscle contraction transforms.

We derive the values for the intracellular and extracellular conductivities needed for bidomain simulations of cardiac electrophysiology using homogenization of partial differential equations. In our model, cardiac myocytes are rectangular prisms and gap junctions appear in a distributed manner as flux boundary conditions for Laplace's equation. Using directly measurable microproperties such as cellular dimensions and end-to-end and side-to-side gap junction coupling strengths, we inexpensively obtain effective conductivities close to those given by simulations with a detailed cyto-architecture (Stinstra et al. in Ann. Biomed. Eng. 33:1743-1751, 2005). This model provides a convenient framework for studying the effect on conductivities of aligned vs. brick-like arrangements of cells and the effect of different distributions of gap junctions along the myocyte membranes.

Of the insects that have been filmed in flight, those that are 1 mm in length or less often clap their wings together at the end of each upstroke and fling them apart at the beginning of each downstroke. This ;clap and fling' motion is thought to augment the lift forces generated during flight. What has not been highlighted in previous work is that very large forces are required to clap the wings together and to fling the wings apart at the low Reynolds numbers relevant to these tiny insects. In this paper, we use the immersed boundary method to simulate clap and fling in rigid and flexible wings. We find that the drag forces generated during fling with rigid wings can be up to 10 times larger than what would be produced without the effects of wing-wing interaction. As the horizontal components of the forces generated during the end of the upstroke and beginning of the downstroke cancel as a result of the motion of the two wings, these forces cannot be used to generate thrust. As a result, clap and fling appears to be rather inefficient for the smallest flying insects. We also add flexibility to the wings and find that the maximum drag force generated during the fling can be reduced by about 50%. In some instances, the net lift forces generated are also improved relative to the rigid wing case.

This article introduces a chemical kinetic model of the transcriptional elongation dynamics of RNA polymerase. The model's novel concept is a look-ahead feature, in which nucleotides bind reversibly to the DNA before being incorporated covalently into the nascent RNA chain. Analytical and computational methods for studying the behavior of the look-ahead model are introduced, and several approaches to parameter estimation are tested on synthetic and also on actual experimental data. Two types of experimental data are considered: 1), the mean velocity of RNA polymerase as a function of the ambient concentrations of the ribonucleoside triphosphates; and 2), the distribution of time intervals between the forward steps of RNA polymerase. By separately fitting the look-ahead model to these two types of data, we obtain estimates of the model parameters. The most difficult parameter to estimate is the width of the look-ahead window. Both types of data suggest a small window size, but the second type does a better job of distinguishing the different window sizes. These latter data rule out a window size of 1, and they strongly suggest a look-ahead window that is approximately four bases in width. Additional experiments to determine the window size are proposed.

We describe and analyze a model for a stochastic pulse-coupled neural network, in which the randomness in the model corresponds to synaptic failure and random external input. We show that the network can exhibit both synchronous and asynchronous behavior, and surprisingly, that there exists a range of parameters for which the network switches spontaneously between synchrony and asynchrony. We analyze the associated mean-field model and show that the switching parameter regime corresponds to a bistability in the mean field, and that the switches themselves correspond to rare events in the stochastic system.

The heart is an electrically controlled fluid pump which operates by mechanical contraction. Whole heart modelling is a computationally daunting task which must incorporate several subsystems: mechanical, electrical, and fluidic. Numerous feedback mechanisms on many levels, and operating at different scales, exist to finely control behaviour. Understanding these interactions is necessary to understand heart operation, as well as pathologies and therapies. A review of the components in such a model is given. The authors then present a framework for their electro-mechano-fluidic whole heart model based on cable methods. The model incorporates atria and ventricles, and has functioning valves with papillary muscles. The effect of altered propagation due to left and right bundle branch block on cardiac output is examined using the cable-based model. Results are compared to clinically observed phenomena. Good agreement was obtained, but tighter coupling of mechanical and electrical events is needed to fully account for behaviour. Cable-based models offer an alternative to continuum models.

We study cardiac action potential propagation under severe reduction in gap junction conductance. We use a mathematical model of cellular electrical activity that takes into account both three-dimensional geometry and ionic concentration effects. Certain anatomical and biophysical parameters are varied to see their impact on cardiac action potential conduction velocity. This study uncovers quantitative features of ephaptic propagation that differ from previous studies based on one-dimensional models. We also identify a mode of cardiac action potential propagation in which the ephaptic and gap-junction-mediated mechanisms alternate. Our study demonstrates the usefulness of this modeling approach for electrophysiological systems especially when detailed membrane geometry plays an important role.

In this paper, we present a one-dimensional model for blood flow in arteries, without assuming an a priori shape for the velocity profile across an artery (Azer, Ph.D. thesis, Courant Institute, New York University, 2006). We combine the one-dimensional equations for conservation of mass and momentum with the Womersley model for the velocity profile in an iterative way. The pressure gradient of the one-dimensional model drives the Womersley equations, and the velocity profiles calculated then feed back into both the friction and nonlinear parts of the one-dimensional model. Besides enabling us to evaluate the friction correctly and also to use the velocity profile to correct the nonlinear terms, having the velocity profile available as output should be useful in a variety of applications. We present flow simulations using both structured trees and pure resistance models for the small arteries, and compare the resulting flow and pressure waves under various friction models. Moreover, we show how to couple the one-dimensional equations with the Taylor diffusion limit (Azer, Int J Heat Mass Transfer 2005;48:2735-40; Taylor, Proc R Soc Lond Ser A 1953;219:186-203) of the convection-diffusion equations to drive the concentration of a solute along an artery in time.

When modulators of neuromuscular function alter the motor neuron spike patterns that elicit muscle contractions, it is predicted that they will also retune correspondingly the connecting processes of the neuromuscular transform. Here we confirm this prediction by analyzing data from the cardiac neuromuscular system of the blue crab. We apply a method that decodes the contraction response to the spike pattern in terms of three elementary building-block functions that completely characterize the neuromuscular transform. This method allows us to dissociate modulator-induced changes in the neuromuscular transform from changes in the spike pattern in the normally operating, essentially unperturbed neuromuscular system.

Individual cells in genetically homogeneous populations have been found to express different numbers of molecules of specific proteins. We investigated the origins of these variations in mammalian cells by counting individual molecules of mRNA produced from a reporter gene that was stably integrated into the cell's genome. We found that there are massive variations in the number of mRNA molecules present in each cell. These variations occur because mRNAs are synthesized in short but intense bursts of transcription beginning when the gene transitions from an inactive to an active state and ending when they transition back to the inactive state. We show that these transitions are intrinsically random and not due to global, extrinsic factors such as the levels of transcriptional activators. Moreover, the gene activation causes burst-like expression of all genes within a wider genomic locus. We further found that bursts are also exhibited in the synthesis of natural genes. The bursts of mRNA expression can be buffered at the protein level by slow protein degradation rates. A stochastic model of gene activation and inactivation was developed to explain the statistical properties of the bursts. The model showed that increasing the level of transcription factors increases the average size of the bursts rather than their frequency. These results demonstrate that gene expression in mammalian cells is subject to large, intrinsically random fluctuations and raise questions about how cells are able to function in the face of such noise., Competing Interests: Competing interests. The authors have declared that no competing interests exist.

The role of protein flexibility in molecular motor function has previously been studied by considering a Brownian ratchet motor that is connected to its cargo by an elastic spring, with the result that the average velocity of the motor/cargo system is increased by reducing the stiffness of the linkage. Here, we extend this investigation to the case of chromosome transport during anaphase A, in which the relevant flexibility is not primarily in the motor/cargo linkage but rather in the cargo itself, i.e., in the chromosome. We model the motor mechanism as an imperfect Brownian ratchet with a built-in opposing load and the chromosome as a collection of discrete segments linked by an elastic energy function that discretizes the potential energy of an elastic rod. Thermal fluctuations are produced in the model by random forces, as in Brownian dynamics. All of the parameters that characterize the chromosome are known or can be estimated from experimental data, as can all but one of the motor parameters, which is adjusted to give the correct transport velocity of normal-length chromosomes. With the parameters so determined, we then reproduce the experimental finding of Nicklas [Nicklas, R. B. (1965) J. Cell Biol. 25, 119-135] that chromosome speed is essentially independent of chromosome length, even though our model contains no "velocity governor." We find instead that this effect is a consequence of chromosome flexibility, as it disappears when stiffer than normal chromosomes are considered.

The kinesin family of motor proteins are involved in a variety of cellular processes that transport materials and generate force. With recent advances in experimental techniques, such as optical tweezers can probe individual molecules, there has been an increasing interest in understanding the mechanisms by which motor proteins convert chemical energy into mechanical work. Here we present a mathematical model for the chemistry and three dimensional mechanics of the kinesin motor protein which captures many of the force dependent features of the motor. For the elasticity of the tether that attaches cargo to the motor we develop a method for deriving the non-linear force-extension relationship from optical trap data. For the kinesin heads, cargo, and microscope stage we formulate a three dimensional Brownian Dynamics model that takes into account excluded volume interactions. To efficiently compute statistics from the model, an algorithm is proposed which uses a two step protocol that separates the simulation of the mechanical features of the model from the chemical kinetics of the model. Using this approach for a bead transported by the motor, the force dependent average velocity and randomness parameter are computed and compared with the experimental data.

Atrial fibrillation (AF) is the most common cardiac arrhythmia. It is characterized by rapid and disorganized electrical activity in the atria. Atrial arrhythmias can be triggered from an ectopic focus, i.e., an abnormal impulse originating in an area other than the sinus node, generating reentrant waves. The regional ionic heterogeneities found in the atria cause a gradual shortening of the action potential duration (APD) with increased distance from the sinoatrial node. It is generally thought that the only electrophysiological consequence of the spatial dispersion of cardiac action potentials (AP) is the enhancement of reentry. This paper investigates the effect of a gradient in APD on arrhythmogenesis via computer simulations. A gradient of ionic properties was introduced into a computationally efficient computer model of the canine atria to produce a smooth distribution of APDs. The window of vulnerability for ectopic beat-induction of reentry was determined for both left atrium (LA) and the right atrium (RA) stimulation, with and without an APD gradient. The shortened windows of vulnerability in the RA, due to the addition of the APD gradient, suggests a protective mechanism against AF. The left atrial window of vulnerability was slightly longer from ionic dispersion.

Circadian (nearly 24-h) clocks are remarkably accurate at timing biological events despite the randomness of their biochemical reactions. Here we examine the causes of their immunity to molecular noise in the context of a detailed stochastic mathematical model of the mammalian circadian clock. This stochastic model is a direct generalization of the deterministic mammalian circadian clock model previously developed. A feature of that model is that it completely specifies all molecular reactions, leaving no ambiguity in the formulation of a stochastic version of the model. With parameters based on experimental data concerning clock protein concentrations within a cell, we find accurate circadian rhythms in our model only when promoter interaction occurs on the time scale of seconds. As the model is scaled up by proportionally increasing the numbers of molecules of all species and the reaction rates with the promoter, the observed variability scales as 1/n(0.5), where n is the number of molecules of any species. Our results show that gene duplication increases robustness by providing more promoters with which the transcription factors of the model can interact. Although PER2 mutants were not rhythmic in the deterministic version of this model, they are rhythmic in the stochastic version.