Click on any of these titles to go directly to the corresponding section: To see publications that appeared during this period: Publications 1998


This report describes work carried out during the 1998 fiscal year in the Mathematical Research Branch, NIDDK. Our work consists of theoretical modeling of systems of biomedical interest, carried out in collaboration with experimental scientists at NIH and elsewhere.

1998 was a year of big transitions for MRB. John Rinzel retired from government service to take a faculty position at New York University, and Arthur Sherman took over as acting chief. Yi-der Chen moved over from NIDDK's Laboratory of Chemical Physical, bringing a new research component in muscle motility and molecular motors. Post-doc Gerda de Vries accepted a faculty position at the University of Alberta, Edmonton, Canada. Sharon Crook successfully defended her thesis and went on to do a post-doc with John Miller and Gwen Jacobs at Montana State University.

MRB hosted a large number of visitors. Robert Miura (UBC, Canada) spent a productive sabbatical year here. Steve Strogatz (Cornell) spent half his sabbatical at MRB. Somdatta Sinha (CCMB, India) came for a short visit and Lee Segel (Weizmann Institute, Israel) spent a month in residence. Hansie Wong (Berkeley) worked as a summer student with Greg Smith.

In eduacational outreach activities, R. Mejia was Visiting Faculty and Colloquium Speaker at the Summer Institute in Mathematics for Undergraduates at the University of Puerto Rico (Humacao Campus) in July, 1998.


A network of oscillatory bursting neurons with excitatory coupling is hypothesized to define the primary kernel for respiratory rhythm generation in the pre-Botzinger complex (pBc) in mammals. Two minimal models of these neurons are proposed, based on different subthreshold ionic current mechanisms for rhythmic bursting. Both models are consistent with many of the dynamic features of electrophysiological recordings from pBc oscillatory bursting neurons in vitro, including voltage-dependent activity modes (silence, bursting, and beating), a voltage-dependent burst frequency which can vary from 0.05 Hz to greater than 1Hz, and a decaying spike frequency during bursting. These results are robust and persist across a wide-range of parameter values for both models. However, the dynamics of Model I are more consistent with experimental data in that the burst duration decreases as the baseline membrane potential is depolarized and the model has a relatively flat membrane potential trajectory during the interburst interval. We propose several experimental tests to demonstrate the validity of either model and to differentiate between the two mechanisms. (RJ Butera, J Rinzel (NYU) and JC Smith (NINDS))

We have proposed models for the ionic basis of oscillatory bursting of respiratory pacemaker neurons in the pBc. We investigated the frequency control and synchronization of these model neurons when coupled by excitatory amino-acid mediated synapses and controlled by convergent synaptic inputs modeled as tonic excitation. Simulations of a heterogeneous population of 50-500 bursting neurons reveal coupling effects similar to those found experimentally in vitro: coupling increases the mean burst duration and decreases the mean burst frequency. Burst synchronization occurred over a wide range of intrinsic frequencies (0.1-1Hz) and even in populations where as few as 10% of the cells were intrinsically bursting. At a population level, both parameter heterogeneity and excitatory coupling synergistically combine to increase the dynamic input range: robust synchronous bursting persisted across a much greater range of parameter space (in terms of mean depolarizing input) than that of a single model cell. This extended dynamic range for the bursting cell population indicates that cellular heterogeneity is functionally advantageous. Our modeled system accounts for the range of intrinsic frequencies and spiking patterns of inspiratory (I) bursting cells found in the pBc complex in neonatal rat brainstem slices in vitro. There is a temporal dispersion in the spiking onset times of neurons in the population, predicted to be due to heterogeneity in intrinsic neuronal properties, with neurons starting to spike before (pre-I), with (I), or after (late-I) the onset of the population burst. (RJ Butera (NIDDK/NINDS), J Rinzel (NYU) and JC Smith (NINDS)).

Multiple Coexisting Periodic Solutions in Modeled Bursting Neurons. Autonomous bursting systems are characterized by periods of repetitive activity punctuated by periods of silence. Recent studies have demonstrated that a model of neuron R15 may possess as many as 8 unique coexisting periodic solutions (multirhythmicity) at a given parameter set. This finding is of novel biological interest, where multirhythmicity has been proposed as a form of short-term memory, and of mathematical interest, since models with more than two coexisting periodic solutions are not frequently encountered. We developed a minimal biophysical model of neuron R15 and a Poincare mapping technique for identifying parameter sets where multirhythmic solutions occur. The origin of multirhythmicity in this class of models is hypothesized and demonstrated in a general 3 variable model, as well as a more complex biophysical model of neuron R15 with 11 state variables. It is possible that other bursting systems with at least two slow variables may exhibit similar forms of multirhythmicity (RJ Butera)

Locomotion Pattern-Generating Properties of a Ring Circuit of Biophysical Neuron Models. A ring circuit of four neuronal oscillators produces patterns corresponding to several quadrupedal gaits, including the walk, bound, and gallop. An analysis using the phase-response curve (PRC) of the component neurons accurately predicted gait modes exhibited by the ring circuit and the phasic relationship between the component neurons. The patterns of network activity are determined by the PRC of the each neuron, which is determined by intrinsic biophysical properties. Altering the PRC of the neurons by a single intrinsic parameter was sufficient to produce gait transitions, with no change in the pattern or strength of connectivity among the neurons. We have also examined the stability of the patterns based on PRC analysis. (RJ Butera, with JW Clark: Rice U., CC Canavier: U. of New Orleans, RO Dror: MIT, and DA Baxter and JH Byrne: U. Texas Med. School at Houston).

Fourier analysis of sinusoidally driven thalamocortical relay neurons and a minimal integrate-and-fire-or-burst model: I have developed a minimal model of the burst and tonic response properties of thalamocortical relay cells observed during sinusoidal current injection in a cat thalamic slice preparation. This ``integrate-and-fire-or-burst'' (IFB) neuron model is constrained using Fourier analysis of both experimental and theoretical firing patterns. The IFB model has many features in agreement with experimental observations such as 1) mixed burst and tonic responses for some stimulus parameters, 2) linearity of tonic responses, and 3) increased phase advance for the nonlinear burst responses. Characterizing the response properties of this minimal model has given insights regarding the stimulus-dependence of burst vs. tonic response mode in TC neurons. (J Rinzel, SM Sherman and L Cox: SUNY Stony Brook, GD Smith)

Spike frequency adaptation in thalamocortical relay neurons and simplified firing rate models: I am developing a firing-rate description of TC cell activity that includes burst responses as well as spike-frequency adaption (SFA) observed during tonic spiking. In this work I have followed XJ Wang's quantitative theory of temporal SFA, originally developed for cortical pyramidal cells [J. Neurophys., 79(3):1549-1566, 1998. The simplicity of such firing-rate neuron models makes them good candidates for large scale network simulations. (J Rinzel, SM Sherman and L Cox: SUNY Stony Brook, GD Smith)

Secretory Physiology

Effect of noise on the emergent collective behaviour of diffusively coupled biological oscillators: We consider a representative model for bursting electrical activity (for example, as observed in pancreatic beta-cells). We choose a parameter regime in which the fast subsystem is monostable so that individual, isolated cells are tonically spiking. Weak diffusive coupling of these cells through gap junctions has been shown to induce stable bursting behaviour. With stronger coupling, the cells revert to tonic spiking. We demonstrate that the addition of noise dramatically increases, via a stochastic resonance-like phenomenon, the coupling range over which bursting is seen. (G de Vries and A Sherman)

Control of action potential firing patterns in GT1 cells: GT1 cells are a cell line derived from hypothalamic, GnRH-secreting neurons. They represent the interface between neuronal and secretory cells and share properties inherent to both groups of cells. These cells are of interest because they represent a key locus in the neuronal control of major endocrine functions, in this case the reproductive axis. This is a new project for the group and we have started by examining the control of spontaneous action potential firing frequency. In particular, our investigations have demonstrated that the cells modulate their firing frequency in response to injection of depolarizing current, i.e. they demonstrate spike frequency adaptation. However, experimental investigations and model analyses demonstrate that typical mechanisms for spike frequency adaptation such as slow, voltage-sensitive ion channels or calcium-activated potassium channels do not regulate spontaneous firing frequency in GT1 cells, so we are attempting to determine what alternative mechanisms may be acting in these cells. (AP LeBeau, A Sherman, and F Van Goor, S Stojilkovic: ERRB/NICHD)


Numerical and analytical results related to localized Ca2+ elevations (Ca2+ ``puffs'' and ``sparks'') viewed as a reaction-diffusion problem: Punctate releases of Ca2+, called Ca2+ sparks, originate at the regular array of t-tubules in cardiac myocytes. I have developed a simple numerical model of Ca2+ spark formation and detection in this cell type. The model involves the numerical solution to a set of reaction-diffusion equations (RDEs) that incorporate Ca2+ release, cytosolic diffusion, and resequestration by SR Ca2+-ATPases, as well as the association and dissociation of Ca2+ with endogenous Ca2+-binding sites and the diffusible indicator dye fluo-3. A realistic feature of these calculations is that the time-dependent spatial profile of fluorescence is convoluted with a point spread function to simulate the contribution of out-of-focus signal. (GD Smith, JE Keizer: UC Davis, MD Stern: NIA, WJ Lederer: UMAB and H Cheng: NIA)

Models of RyR-mediated Ca2+-induced Ca2+ release and the spark-to-wave transition in cardiac myocytes: During Ca2+ overload, Ca2+ sparks serve as sites for the initiation of Ca2+ waves in myocytes. I have carried out computer simulations of spark-induced waves to explore the influence of the regular array of release sites on their propagation. These computer simulations with a linear array of Ca2+ release sites give a wave speed proportional to the Ca2+ diffusion constant rather than its square root, as is true for classical traveling waves in an excitable medium. (JE Keizer: UC Davis, GD Smith, S Ponce-Dawson and J Pearson: Los Alamos)

Mathematical analysis of reaction-diffusion equations (RDEs) associated with the buffered diffusion of Ca2+: Analysis of the RDEs describing the buffered diffusion of Ca2+ in cytosol reveals that two dimensionless parameters, one dependent on magnitude of Ca2+ release, determine whether the cytosol is in the ``excess buffer'' or ``rapid buffer'' regime. We are performing regular and singular perturbation analysis leading to higher order approximations in both limits. (L Dai, R Miura, GD Smith, A Sherman.)

Biological Free Energy Transduction

Theoretical formalisms for motor motility: Biological motors, such as kinesins and dyneins, are proteins that can catalyze the hydrolysis of ATP and use the released energy to move a load (such as a vesicle in axonal transports) on a microtubule filament. In vitro, both the biochemical catalysis and the mechanical movement of single motors have been studied separately. However, exactly how the two processes are coupled is still not clear. The reason for this is that formalisms relating the catalytic cycle of ATP hydrolysis and the movement of the motor on a microtubule are still lacking. This project is devoted to solving this problem. We first consider the well-studied in vitro assay in which a large plastic bead is attached elastically to a single motor molecule and the movement of the bead is studied in a microscope. Last year, a theoretical formalism was obtained for the case that the motor was one-headed. The formalism can be easily used to calculate the velocity of the bead in the assay when the parameters of the biochemical cycle are given. This year, we apply the formalism to a simple three-state catalysis model to examine how the velocity of the movement of the bead in this assay is affected by the viscosity of the medium, the stiffness of the elastic element connecting the bead and the motor, etc. We have found that the model can easily generate the mechanical properties of kinesin motors measured experimentally. (Chen).

Regulation of actin-myosin interaction by caldesmon: Muscle contraction is generated by binding and unbinding of myosin molecules to actin while catalyzing ATP hydrolysis. Caldesmon, a protein exists in smooth muscles, can bind to actin and reduce the ATP hydrolysis rate. How caldesmon regulates the actin-activated ATP hydrolysis of myosins is still not known. One theory is that caldesmon acts as a pure inhibitor for myosin binding. Last year, we showed that this might be the case based on kinetic measurements and Monte Carlo simulations. In that study, it was assumed that the number of sites on actin covered by a bound caldesmon was seven and that the actin can exist only in one state. These assumptions are considered to be reasonable previously, but not recently. Thus, we carried out a series of systematic model calculations (Monte Carlo simulations) to study the effect of these two assumptions. We found that the simultaneous binding of myosin and caldesmon to actin may not be pure competitive as found previously. (Yan, Chen, and J. M. Chalovich at University of East Carolina).

Fluctuation-induced biased Brownian motion: Brownian particles can be made to execute biased movement on a one-dimensional periodic ratchet-like potential if the potential is turned on-and-off randomly or regularly and each ratchet is asymmetric. In general, the direction of the biased movement is independent of the frequency of the fluctuation of the potential with regular ratchets. In this study, we found that frequency-modulated direction reversal can be found in this system if the ratchet is distorted by kinking one arm vertically. The results should be useful in designing the apparatus for separating biomolecules based on their sizes and charges. (Chen, Yan, and R. Miura (University of Vancouver, Canada))


The hypothesis used is that the vasomotion pattern depends on a metabolic parameter (Berne) that results from the tissue oxygen availability and the tissue metabolic demand level, and that influences the terminal arteriole dynamics. The model incorporates the terminal arteriole and capillary bed structures, oxygen transport, tissue diffusion, tissue energy metabolism (ATP lysis, mitochondria oxidative phosphorylation, and ATP reconversion), arteriolar wall ionic transports, its muscle-elasticity dynamics, and blood flow. The model shows that as the metabolic demands increases the blood flow increases. At moderate and high metabolic demand levels this adjustment results from a progressive lengthening of the vasorelaxation phase and shortening of the vasocontraction phase. On the other hand, for basal and low metabolic demand levels the accompanying vasomotion patterns are complex and at some instances chaotic, yet the time average blood flow increases pari pasu with the metabolic demand level. The latter may correspond wit h the non regular vasomotion (Bassingwaithe) adjustments observed experimentally (Jose M. Gonzalez-Fernandez, Arthur Sherman, MRB, NIDDK; Bard Ermentrout, U. Pittsburgh).

Structural Biology

Integral membrane proteins comprise receptors and transporters in epithelial cells that form the mammalian kidney. Neural net models that are trained using well-characterized proteins (e.g., G-protein coupled receptors) have been developed to identify the location of protein segments relative to the membrane or cytosole. Several proteins of importance in renal function have been studied, including aquaporins (AQP1, AQP2, AQP3, AQP4), the bumetanide-sensitive Na-K-2Cl co-transporter (BSC-1), and urea transporters (UT1, UT2). (R Mejia and M Knepper: NHLBI)

Renal Physiology

Urine concentration mechanism. Three major clinical states that manifest water and NaCl retention in their severe forms are congestive heart failure, cirrhosis, and nephrotic syndrome. Mathematical modeling studies are being conducted to determine mechanisms of water retention that might account for development of hyponatremia. (R Mejia, P Fernandez-Llama and M Knepper: NHLBI)

Concentration mechanism in the renal inner medulla. The mechanism by which the mammalian kidney concentrates solutes in the inner medulla has yet to be established, in particular, because active solute transport is considered unlikely due to limited blood flow. We have developed models to test several hypotheses. These have shown that an elastic epithelial cell with two bathing solutions and distinct apical and basolateral transport properties will not concentrate fluid passively under peristaltic stress. A model of the inner medulla that includes both cellular and tubular compartments shows modest gradients due to peristalsis and larger gradients obtained with the cellular compartment. This is due to osmotic driving forces that develop in the cell. (R Mejia, M Knepper: NHLBI,

Immunomorphometry. The anatomy and histology of the kidney have been studied for decades, most recently in an effort to understand the "passive concentrating mechanism". We use immunofluorescent immunolabeling of renal tissue sections to identify both structure and function. Transporters being used to identify and quantify structure include AQP1 (located in descending Henle's limb (DL)), AQP2 (located in collecting duct (CD)), UT (located in both DL and CD); a Cl- channel protein (located in ascending Henle's limb); and von Willebrand factor located in the vasa recta. (R Mejia, J Wade: UMD,, M Knepper: NHLBI)

Parameters for experimentalists and modelers. A database of renal parameters with author, title, abstract, address and source continues to be compiled. Physiologic and geometric parameters as well as conditions of measurement are included. This database, which is in extended MEDLINE format, was formerly made available to investigators via a gopher server. This server was discontinued early in 1998, and alternatives for availability of the database on the web are being considered; currently parameter references are available upon request to (R Mejia and M Knepper: NHLBI)

Cell Energetics

We have previously described the effect of changing transport work on the concentration profiles of high energy phosphate compounds within single cells using a reaction-diffusion model. A model with cylindrical geometry (pdes in two space dimensions and time) is used to investigate the biological hypothesis that ATP concentration is not limiting and that there can be significant ADP concentration gradients within the cell. We presently test the hypothesis that there are local changes near the membrane that modulate channel function without significant changes in global cystolic concentration, and we do this by incorporation of vesicular compartments into the model. (R Mejia and R Lynch: U Arizona,


Measurement of agonist and antagonist ligand binding parameters at the human parathyroid hormone type 1 receptor: The aim of this study was to develop radioligand binding assays for the PTH-1 receptor in order to measure agonist and antagonist binding parameters. Agonist and antagonist binding were compared to identify characteristics that are specific for ligands that activate the receptor. Ligand binding was compared in the presence and absence of 5'-O-(3-thiotriphosphate)(GTPgammaS) to evaluate the effects of G-protein coupling on receptor-ligand interaction. Data obtained in the presence of GTPgammaS were used to address ligand binding properties at the G-protein-uncoupled receptor. The data were used to evaluate potential models of receptor-ligand interaction. The most plausible explanation assumes two binding sites on the receptor that recognize two corresponding sites on agonist ligands. (SRJ Hoare:NIMH, G de Vries, and TB Usdin:NIMH)

Please send your comments/suggestions about this page to:

To go back to MRB's Home Page click here