A biochemical description of postsynaptic plasticity-with timescales ranging from milliseconds to seconds.

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.

Strong intracellular signal inactivation produces sharper and more robust signaling from cell membrane to nucleus.

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.

Mathematical modeling of the impact of recirculation on exchange kinetics in tandem extracorporeal membrane oxygenation and therapeutic plasma exchange.

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.

Spontaneous oscillation and fluid-structure interaction of cilia.

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.

Modeling and high-throughput experimental data uncover the mechanisms underlying Fshb gene sensitivity to gonadotropin-releasing hormone pulse frequency.

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.

Image-based model of the spectrin cytoskeleton for red blood cell simulation.

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.

Improved signaling as a result of randomness in synaptic vesicle release.

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.

A Mathematical Model of Granule Cell Generation During Mouse Cerebellum Development.

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 influence of volume exclusion by chromatin on the time required to find specific DNA binding sites by diffusion.

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.

Computer simulation of voltage sensitive calcium ion channels in a dendritic spine.

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.

The influence of spatial variation in chromatin density determined by X-ray tomograms on the time to find DNA binding sites.

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.

Quantitative theory of telomere length regulation and cellular senescence.

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.

Computer simulation of surgical interventions for the treatment of refractory pulmonary hypertension.

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.

Simulating Cardiac Fluid Dynamics in the Human Heart.

Cardiac fluid dynamics fundamentally involves interactions between complex blood flows and the structural deformations of the muscular heart walls and the thin, flexible 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. However, existing computational models that account for interactions among the blood, the actively contracting myocardium, and the cardiac valves are limited in their abilities to predict valve performance, resolve fine-scale flow features, or use realistic descriptions of tissue biomechanics. Here we introduce and benchmark a comprehensive mathematical model of cardiac fluid dynamics 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 fluid-structure interaction 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 predicts fine-scale flow features. None of these outputs are prescribed; instead, they emerge from interactions within our comprehensive description of cardiac physiology. Such models can serve as tools for predicting the impacts of medical devices or clinical interventions. They also can serve as 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.

Synchrony and asynchrony for neuronal dynamics defined on complex networks.

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 (Erdos-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.

The impact of universal recycling on the evolution of economic diversity.

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.

Optimal Fenestration of the Fontan Circulation.

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.

Ephaptic conduction in a cardiac strand model with 3D electrodiffusion.

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.

Models for plasma kinetics during simultaneous therapeutic plasma exchange and extracorporeal membrane oxygenation.

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.