Second Young Reseachers Workshop in Mathematical Biology: Abstracts and Lecture Materials
Plenary Speakers
Animals, including humans, use interaural time differences (ITDs) that arise from of different sound path lengths to the two ears, as a cue of horizontal sound source location. The nature of the neural code for ITD is still controversial. Current models advocate either a map-like place code of ITD along an array of neurons, consistent with a large body of data in the barn owl, or a rate-based population code, consistent with data from small mammals. Recently, it was proposed that these different codes reflect an optimal coding strategy that depends on head size and sound frequency. The chicken makes an excellent test case because its physical prerequisites are similar to small mammals, yet it shares a more recent common ancestry with the owl. We show here that, like in the barn owl, the brainstem nucleus laminaris in mature chickens displayed the major features of a place code of ITD. The physiological range of ITDs was systematically represented in the maximal responses of neurons along each isofrequency band. This is in contrast to the predictions from optimal coding theory and thus re-opens the question as to what determines the neural coding strategies for ITDs, including which code might be implemented by the human brain.
I will describe some recent work in our group on the dynamics of the actin cytoskeleton in relation to the movement of a motile cell. First, I will describe work (joint with Adriana Dawes, Eric Cytrynbaum, and Bard Ermentrout) on a simple 1D spatial model of a cell. We show how the branching of actin filaments and the forces they exert on the cell membrane account for the protrusion velocity and characteristic actin density profiles. (This work is partly analytical and partly numerical.) We use this model to understand how branching rates and other biochemical parameters control cell speed by studying the relevant travelling wave solutions.
I will also describe efforts (joint with AFM Maree, Alexandra Jilkine, Adriana Dawes and Veronica Grieneisen) at assembling a more detailed 2D spatial model of a crawling cell, in which we take into account the regulatory role of a set of signalling proteins (Cdc42, Rac, Rho). We show how the interplay between these and the actin cytoskeleton accounts for the ability of the cell to self-organize, polarize, maintain a stable shape and speed, and respond to new external signals.
What aspects of a stimulus cause a neuron to fire? How do stimuli affect the time of spikes? In this talk, I will discuss what we can learn about neuronal firing patterns by regarding neurons as nonlinear oscillators. The spike-triggered average or reverse correlation method is a common approach for determining what kinds of stimuli make a neuron fire. The poststimulus time histogram is another experimental measurement for describing the affect of a stimulus on the firing pattern of a neuron. The latter can be related to the former by using some optimality arguments. Both of these curves should be affected by the membrane properties of the individual neuron of interest. Since this is a huge-dimensional space, we will focus on one property of neurons which has been shown to be tightly coupled to neuronal dynamics: the phase resetting curve (PRC). The PRC describes the shift in the timing of a spike due to a brief stimulus as a function of the time since the last spike. We show that under certain circumstances there is a 1:1 mapping between the STA, the PSTH, and the PRC. Thus, we connect internal dynamics of neurons with their preferred stimuli and their population responses. This work is joint with Boris Gutkin, Alex Reyes, Nathan Urban, Roberto Galan, and Nicolas Fourcaud.
This lecture will be devoted to a discussion of some of the basic problems in modeling actin dynamics. Actin polymerization and network formation are key processes in cell motility. Numerous actin binding proteins controlling the dynamic properties of actin networks have been studied and models such as the dendritic nucleation scheme have been proposed for the functional integration of at least a minimal set of such regulatory proteins. However, a complete understanding of actin network dynamics is still lacking. Even at the actin-filament level, the dynamics of the distribution of filament lengths and nucleotide profiles are still not fully understood. We will describe recent work on the evolution of the distribution of filament lengths and nucleotide profiles of actin filaments, both from a deterministic and a stochastic viewpoint. If time permits we will discuss work aimed at integrating microscopic models of actin dynamics into cell-level descriptions of motility.
The main function of the circulatory system is to transport and exchange substances throughout the body. Delivery of oxygen is a particularly demanding function, because oxygen is relatively insoluble in water. Within blood vessels, oxygen is carried convectively by hemoglobin molecules within red blood cells. Oxygen exchange with tissue occurs by diffusion in the microcirculation, an extensive branching network of microscopic vessels that brings blood close to all oxygen-consuming tissues. The microcirculation regulates blood flow according to changing local demands over short and long time scales. Mathematical models can be used to gain insight into these processes. Models will be described for the mechanics of blood flow in capillaries, for oxygen exchange between blood and tissues and for structural adaptation of blood vessels. Applications to disease states including cancer will be discussed.
The first generation of models for electrical activity in pancreatic beta-cells focused on ionic mechanisms. Negative feedback by calcium, directly onto calcium-activated potassium channels and indirectly onto ATP-sensitive potassium channels and sodium pumps, is the main type of mechanism considered in current models. Such models do a good job of accounting for the oscillations on a wide range of time scales, ranging from 10 seconds to about 2 minutes. However, even slower oscillations, with periods of 4 or even 10 minutes are often observed, and these often appear with the faster oscillations layered on top. This suggests that there is an additional mechanism for oscillations, which we have proposed is based on oscillations of glycolysis. We will discuss how two relatively simple oscillators, which can be off, oscillating, or tonically on, depending on stimulation level, can be combined to account for the great diversity of observed patterns. We will also consider the impact of metabolic oscillations on synchronization of beta-cells within the islet of Langerhans. Diffusion of glycolytic metabolites provides an important mechanism for secretion, but can also lead to oscillator death and a source of bistability.
Gliomas account for over half of all primary brain tumors and have been studied extensively for decades. Even with increasingly sophisticated medical imaging technologies, gliomas remain uniformly fatal lesions. A significant gap remains between the goal of designing effective therapy and the present understanding of the dynamics of glioma progression. It has become increasingly clear that, along with the proliferative potential of these neoplasms, it is the subclinically diffuse invasion of gliomas that most contributes to their resistance to treatment. That is, the inevitable recurrence of these tumors is the result of diffusely invaded but practically invisible tumor cells peripheral to the abnormal signal on medical imaging and to the limits of surgical, radiological and chemical treatments.
In this presentation, I will demonstrate how quantitative modeling can not only shed light on the spatio-temporal growth of gliomas but also can have specific clinical application in real patients. Integration of our quantitative model with the T1-weighted and T2-weighted magnetic resonance (MR) imaging characteristics of gliomas can provide estimates of the extent of invasion of glioma cells peripheral to the imaging abnormality. Additionally, further insight can be gained from parametric mapping of kinetic model parameters derived from positron emission tomography (PET) with novel tracers. In summary, although current imaging techniques remain woefully inadequate in accurately resolving the true extent of gliomas, quantitative modeling provides a new approach for the dynamic assessment of real patients and helps direct the way to novel therapeutic approaches.
Poster Presentations
Nuclear Magnetic Resonance (NMR) is a spectroscopic technique that involves the study of molecular structures via the interaction of radio-frequency electromagnetic radiation with the nuclei immersed in a strong magnetic field. NMR experiments provide typically two kinds of structural constraints; distance constraints using solution state NMR and orientational constraints using solid state NMR (ssNMR). Membrane protein structures are extemely hard to determine. ssNMR has been found to be effective in solving transmembrane protein structures. The mathematical problems that arise while solving for the three dimensional atomic structure of a membrane protein from ssNMR data are discussed here. Tools from linear algebra and differential geometry are used.
A new variational region based model for a simultaneous image segmentation and registration is proposed. The purpose of the model is to segment and register novel images simultaneously using a modified Mumford-Shah technique and region intensity values. The segmentation is obtained by minimizing a modified Mumford-Shah model. A global rigid registration is assisted by the segmentation information and region intensity values. In addition, the model can also be applied to the case of non-rigid registration. The numerical experiments of the proposed model are tested against simulated normal noisy human-brain MRI images. The preliminary experimental results show the effectiveness of the model in detecting the boundaries of the given objects and registering novel images simultaneously.
Evidence from microscopy suggests that the spatial organization and heterogeneity of receptors and other signaling proteins is a key feature of many signaling pathways. Yet, the vast majority of mathematical models for cell signaling (ODE-based reaction network models) are based on the assumption that each biochemical species is uniformly distributed throughout the cell. The widespread use of ODEs is in part due to the increased computational cost of spatiotemporal modeling and in part due to the lack of detailed knowledge of the model parameters associated with the passive and active transport mechanisms of the cell. While the spatiotemporal modeling is more challenging than well-mixed modeling, spatiotemporal modeling can provide new insights into microscopic events leading to an improved understanding of signal transduction processes that are important at the molecular level. The work presented here couples kinetic Monte Carlo simulations (used to describe the spatial and temporal dynamics of the ErbB receptor on the cell membrane) with reaction-advection-diffusion systems of PDEs (used to describe downstream ErbB intracellular signaling). The results clearly show that spatial heterogeneity is important for the control and regulation of the ErbB network.
Bifurcation theory and perturbation theory can be combined with a knowledge of the underlying circuitry of the visual cortex to produce an elegant story explaining the phenomenon of visual hallucinations. A key insight is the application of an important set of ideas concerning spontaneous pattern formation introduced by Turing in 1952. The basic mechanism is a diffusion driven linear instability favoring a particular wavelength that determines the size of the ensuing stripe or spot periodicity of the emerging spatial pattern. Competition between short range excitation and longer range inhibition in the connectivity profile of cortical neurons provides the difference in diffusion length scales necessary for the Turing mechanism to occur and has been proven by Ermentrout and Cowan to be sufficient to explain the generation of a subset of geometric hallucinations reported. Incorporating further details of the cortical circuitry, namely that neurons are also weakly co! nnected to other neurons sharing a particular stimulus orientation or spatial frequency preference at very long ranges and the resulting shift-twist symmetry of the neuronal connectivity, improves the story. We expand this approach in order to be able to include the tuned responses of cortical neurons to additional visual stimulus features such as motion, color and disparity. We apply a study of nonlinear dynamics similar to the analysis of wave propagation in a crystal lattice to demonstrate how the spatial pattern formed through the Turing instability can be pinned to the geometric layout of various feature preferences. The perturbation analysis is analogous to solving the Schrödinger equation in a weak periodic potential. Competition between the local isotropic connections which produce patterns of activity via the Turing mechanism and the weaker patchy lateral connections that depend on a neuron's particular set of feature preferences create long wavelength affects analo! gous to commensurate-incommensurate transitions found in fluid systems under a spatially periodic driving force. In this way we hope to better understand how the intrinsic architecture of the visual cortex can generate patterns of activity that underlie visual hallucinations.
The challenge for drug design is to create molecules with optimal function that also partition into appropriate in vivo compartment(s). Pharmacokinetics can be improved by ancillary molecules, such as cyclodextrins, tha increase the effective concentrations of hydrophyic drugs in the blood by encapsulating the drug. However, many drug targets are located inside cells and encapsulating molecules have not been developed that would increase the effective concentration of drugs inside cells. Here we propose that RNA aptamers might preform the same functions inside cells, but with higher specificity, as do cyclodextrins in the body fluids.
Existing ecological theory suggests that environmental noise should increase the amount of variation observed in the size of populations across time. In contrast, two recent empirical studies suggest that populations may exhibit less variation in a more variable environment. Here I present theory illustrating 4 mechanisms that could explain this result. Environmental variation can reduce the amplitude of cycles due to 1) non-linear averaging, 2) counter-acting other environmental variation, 3) destructive interference, or 4) compensatory dynamics at a tropic level. In real systems, these mechanisms may act in separately or in concert.
Plants have two main modes of pollen dispersal; via wind or by animal. To date, both of these mechanisms have been examined using the same diffusive models. In order to better understand the subtleties of these dispersal modes, a cellular automata model is used to predict pollen range, genetic diversity, and genetic differentiation. Among the main controlling factors in the model are propensity on the pollen movement, the probability of pollination, and the density of individuals of the species.
We present an automated spectral interpretation algorithm that determines peptide sequence identity directly from MS tandem mass spectra (de novo sequencing). This technology and its ability to perform de novo sequence identification, independent of a genome sequence or annotation, will be of critical importance for localization of single nucleotide polymorphisms, identification and characterization of translational editing, and the detection of genes where either genomic sequences are not available or their annotation is incorrect.
Data discretization, also known as binning, is a frequently used technique in computer science, statistics, and their applications to biological data analysis. A new discretization method is presented which was motivated by a new reverse engineering approach to modeling the biochemical network governing the cellular response of yeast under oxidative stress. While the method can be used for data clustering, it was developed to answer the need for a discretization algorithm that is specifically designed to handle typical systems biology data: a large number of short multivariate time series of heterogeneous data, such as transcript, protein, and metabolite concentration measurements. To such data sets, statistical discretization methods are hardly applicable due to the prohibitive cost of obtaining sample sets of sufficient size. The method transforms real-valued data into a finite number of discrete values. Novel aspects are the incorporation of an information-theoretic criterion and a criterion to determine the optimal number of values. As several modeling techniques for biochemical networks employ discrete variable states, the method needs to preserve the dynamic features of the time series as well as be robust to noise in the experimental data. An example demonstrating the method's qualities will be shown.
We consider neuronal models where the transition from resting behavior to spiking behavior takes place via a Hopf Bifurcation. These models can exhibit resonance, where a certain frequency of input elicits the largest response in the cell. We study the case where the synaptic input to such a cell undergoes both facilitation and depression. The result is a band pass filter where the synapse is strongest at a certain driving frequency. This is known as synaptic resonance. We look for driving frequencies for the synapse that result in the maximum output in the resonant cell.
We consider the directed evolution of a population after an intervention that has significantly altered the underlying fitness landscape. We model the space of genotypes as a distributive lattice; the fitness landscape is a real-valued function on that lattice. The risk of escape from intervention, i.e., the probability that the population develops an escape mutant before extinction, is encoded in the risk polynomial. Tools from algebraic combinatorics are applied to compute the risk polynomial in terms of the fitness landscape. In an application to the development of drug resistance in HIV, we study the risk of viral escape from treatment with the protease inhibitors ritonavir and indinavir.
The introduction of delays into ordinary differential equation models of populations complicates the dynamics considerably. Often, globally attractive equilibria become unstable, and periodic solutions emerge. I study this scenario in one and two dimensional population and predator-prey models.
Prediction of protein residue-residue contacts from sequence information is an important and difficult problem in structural bioinformatics. Most methods use correlated mutation analysis to detect such contacts. We developed a new approach (P2PConPred) for intra-protein contact prediction, which is oriented for detecting direct physical contacts. Our method uses a novel pair-to-pair substitution matrix (P2PMAT) derived from accurate protein multiple sequence alignments with available representative structures. The P2PMAT matrix integrates the probabilities of contacting and non-contacting protein sites. Incorporating evolutionary (sequence) conservation of residues with information regarding correlated substitutions is natural and effective in our P2PMAT matrix. Our P2PConPred method is most sensitive for contact prediction in the protein cores. Core residues are effectively identified from sequence information alone by means of a predicted surface accessibility of! proteins. Our method improves protein core contact prediction by 1.25 and 1.6 fold over the contact prediction method of Gobel et. al. and that of Singer et. al., respectively. Combining our approach with other approaches for calculating correlated mutations is expected to be beneficial. The basic approach we developed can also be naturally applied to protein structure prediction and de-novo drug design.
The Anchor Cell (AC) and the Ventral Uterine Precursor Cell (VU) participate in a critical interaction during the development of the C.elegans hermaphrodite gonad. We propose a new mechanism for this interaction and show, via analysis and computation of our mathematical model, that this mechanism resolves outstanding questions regarding the AC/VU fate decision.
Cardiac mechanics can be modeled as the dynamic interaction of a viscous incompressible fluid (the blood) and a (visco-)elastic structure (the muscular walls and the valves of the heart). The immersed boundary (IB) method is a mathematical formulation and numerical approach to such problems. In this presentation, we describe an adaptive version of the IB method introduced in [1,2]. The adaptive scheme employs the same hierarchical structured grid approach (but a different numerical scheme) as the two-dimensional adaptive IB method of Roma et al. [3,4] and is based on a formally second order accurate (i.e., second order accurate for problems with sufficiently smooth solutions) version of the IB method that we have recently described [5]. We present results obtained from the application of this adaptive methodology to the three-dimensional simulation of blood flow in the heart. The results obtained by the adaptive method show good qualitative agreement with simu! lation results obtained by earlier non-adaptive versions of the method, but the flow in the vicinity of the model heart valves indicates that the new methodology provides dramatically enhanced boundary layer resolution. Differences are also observed in the details of the flow about the mitral valve leaflets.
- B.E. Griffith, Simulating the blood-muscle-valve mechanics of the heart by an adaptive and parallel version of the immersed boundary method, Ph.D. thesis, Courant Institute of Mathematical Sciences, New York University (2005).
- B.E. Griffith, R.D. Hornung, D.M. McQueen, C.S. Peskin, An adaptive, formally second order accurate version of the immersed boundary method, Submitted.
- A.M. Roma, A multilevel self adaptive version of the immersed boundary method, Ph.D. thesis, Courant Institute of Mathematical Sciences, New York University (1996).
- A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (2) (1999) 509--534.
- B.E. Griffith, C.S. Peskin, On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems, J. Comput. Phys. 208 (1) (2005) 75--105.
Stoichiometric Network Theory is a constriants-based, optimization approach for quantitative analysis of the phenotypes of large-scale biochemical networks while avoiding the use of detailed kinetics. This approach uses the reaction stoichiometric matrix in conjunction with constraints provided by flux balance and energy balance to guarantee mass conserved and thermodynamically allowable predictions, respectively. However, to date, the flux and energy balance constraints have not been implemented simultaneously because optimization under the combined constraints is nonlinear. I will present a sequential quadratic programming algorithm that solves the nonlinear optimization problem. A simple example and the system of fermentation in Saccharomyces cerevisiae are used to illustrate the new method. The algorithm allows the use of nonlinear objective functions. As a result, a novel optimization with respect to the heat dissipation rate of a system can be suggested. The importance of incorporating interactions between a model network and its surroundings is also emphasized.
Software which is readily available to the neuroscience community can be used to reconstruct cortical surfaces from magnetic resonance imaging (MRI) data. Topologically correct cortical surfaces representing a white matter (WM) surface (which occurs at the white matter/gray matter interface) and a gray matter (GM) surface (which occurs at the gray matter/cerebrospinal fluid interface) can be created. Using these surfaces, we present an approach for producing a cortical mantle volume representing the gray matter. We then use this cortical mantle volume to restrict the analysis of functional MRI (fMRI) data to the gray matter. These results are then compared to a similar analysis using the entire masked volume.
The medial entorhinal cortex (mEC) participates in theta (4-12Hz), beta (15-25Hz), and gamma (30-80Hz) frequency rhythms, which are thought to play an important role in the formation of neuronal ensembles. In vivo studies have shown theta-nested gamma oscillations in the superficial layers of the EC during exploratory behavior and REM sleep, while in-vitro studies have shown theta modulation of gamma activity during kainate application. In vivo, the theta component is thought to arise mainly due to inputs from the medial septum. Recent findings from slice preparations have shown several rhythmic activity patterns in the superficial layers of the entorhinal cortex during the presence and absence of kainate and NMDAR antagonists. Under kainate, NMDAR block causes a decrease in gamma activity, while in the absence of kainate gamma activity increases. We investigate the mechanisms for these surprising experimental results with a model consisting of populations of pyra! midal cells, fast-spiking interneurons, theta-producing interneurons, and stellate cells. In addition, we investigate rhythms generated by the addition of septal input to the model.
In this study, we combine insights from human macroscopic experimental measures and computational neural network modeling to test the hypothesis that perception correlates with cortical activity measured in both the time and frequency domain, and that this activity is mediated by specific cellular level neuronal events. This is accomplished using a two-fold approach. First, we experimentally probe the influence of perception on human cortical rhythms using a somatosensory tactile detection task. Specifically, we will use techniques recently developed at the Arthinoula A. Martinos Center for Biomedical Imaging to measure magnetoencephalography (MEG) signals during presentation of a threshold level somatosensory stimulus. We analyze the signals generated in the primary somatosensory cortex (SI). In the time domain, we will measure amplitudes and latencies of evoked responses, and in the frequency domain will measure spectral power and phase-locking. We will compare ! these measures during perceived and unperceived trials. Second, we use biophysically based neural network modeling to test if changes in the level of " top down" modulation create a systematic biophysical link between changes in time and frequency domain activity, which are consistent with the experimental results. This approach entails the development of a model of a laminated cortical column(s) that reproduces the oscillatory current dipoles that are measured extracranially with MEG. This two-fold approach may lead to a better understanding of the macroscopic and cellular mechanisms of perception.
Cell motility is essential to cancer invasion. Cells use a variety of directed migration mechanisms in response to chemokinesis (chemotaxis), rigidity (mechanotaxis), or extracellular matrix gradients (haptotaxis). Many cancer studies focused on chemotaxis as the main event for a cell to invade, however, little is known on the role of haptotaxis in cancer invasion. Our goal is to produce a quantitative understanding of the haptotatic component of cancer invasion by using a mathematical model. To this end, our first step is to implement original techniques for measuring physical parameters of haptotaxis under controlled conditions. We modified a methodology developed at Vanderbilt to study chemotaxis (Walker et al., 2005). Instead of creating a chemokine gradient, we used a gradient mixer to lay down a gradient of fluorescently-labeled extracellular matrix protein. Then, cells were introduced into the device and allowed to adhere and spread. Haptotaxis was followed ! by videomicroscopy. Fluorescent pictures were acquired to assess the slope of the gradient. Brightfield movies were analyzed to calculate parameters such as cell speed and trajectory.
Our results show that we could efficiently and reproducibly prepare matrix gradients using our device, which is suitable for live-cell imaging. Moreover, using different matrix proteins, we were able to demonstrate the usefulness of this methodology to produce reliable and reproducible migration substrates. We also demonstrated that cells respond differently to matrix gradients with distinct slopes and/or concentration ranges. The methodology we present here is efficient in quantitatively measuring haptotaxis on various matrix gradients. In the future, we will develop gradients of different matrix or even controlled heterogeneous matrix pattern to correlate cancer cell invasion with haptotaxis abilities of invasive cells.
Excitatory-inhibitory networks arise in many neuronal systems. Examples include models for thalamic sleep rhythms and parkinsonian tremor. Such networks have been shown to exhibit a rich structure of firing patterns, including synchronous activity, irregular and chaotic dynamics and propagating wave-like behavior. Computational and analytical methods have been employed extensively to understand those patterns in networks with one spatial dimension. However there has been very little work devoted to the numerical or analytical investigation of higher dimensional networks. We consider two-dimensional sheets of synaptically coupled excitatory and inhibitory neurons and explore the types of additional patterns that emerge. The models consist of large systems of nonlinear differential equations and represent the interactions between two neural populations: the subthalamic nucleus and the globus pallidus. The membrane potential in those models exhibits bursting patterns and thus reveals several time scales. Using a dynamical systems approach we analyze the mechanisms underlying such bursts and then reduce this complex high-dimensional model to a simpler yet biophysically meaningful system.
In the second part of this project, We use the reduced model to derive conditions on network parameters for the existence of various propagating patterns. We compute the functional dependence of the velocity on parameters controlling the inhibitory synaptic input and explain the failure of propagation that occurs for a certain parameter range.
In the current study, we derive a mathematical model describing fluorescence recovery after photobleaching (FRAP) of the Golgi complex, the endoplasmic reticulum (ER) or cytosol in terms of interactions between different pools of depalmitoylation mutant GFP-HRas proteins (DMGRs).
Governing equations for each pool of GFP-Ras proteins are derived as a system of partial differential equations (PDEs) and all the kinetic constants are obtained by Repetitive Fluorescence Recovery after Photobleaching (ReFRAP), a new way to study binding kinetics in vivo. Using the model constructed, we show that (1) fluorescence intensity indicates not only the concentration of fluorescent molecules but also the density of membrane to which fluorescent molecules bind, (2) both kinetics of bleached DMGRs and active fluorescent DMGRs independent if partitioning type binding kinetics is assumed, (3) DMGRs are maintained at dynamic equilibrium by diffusion and binding kinetics, but in similar time scales (4) the bindings of DMGRs on both ER and Golgi membranes follow partitioning kinetics rather than ligand-receptor type binding kinetics, (5) GFP-Ras protein binding kinetics on the ER and the Golgi membrane are distinct with different binding rate constants and (6) ReFRAP can ! be used to determine binding rate constants in both partitioning and ligand-receptor types bindings. And finally (7) it will be also shown that the model has unique steady state which is globally asymptotically stable.
When paced at a rapid rate, cardiac tissue typically undergoes a period-doubling bifurcation to a 2:2 alternans rhythm, where successive action potentials alternate between having long and short duration. During discordant alternans different spatial regions alternate out of phase. Alternans can induce large repolarization gradients across the heart, providing a substrate for unidirectional block and ensuing cardiac arrhythmias; therefore efforts into understanding how alternans arises are important. Structural barriers to wave propagation in cardiac tissue are associated with a decreased threshold for alternans both experimentally and clinically. Using computer simulations, we investigated the effects of a structural barrier on the onset of spatially concordant and discordant alternans. We used a two-dimensional tissue geometry with anisotropy, as well as gradients in selected potassium channel densities to mimic known apex-base gradients. In ionically homogeneous! tissue, discordant alternans arises at intermediate pacing rates due to a source-sink mismatch behind the barrier. The structural barrier also facilitates discordant alternans in ionically heterogeneous tissue. We determined the mechanism of discordant alternans to be a barrier-induced decoupling of tissue with different restitution properties. Our results demonstrate a causal relationship between the presence of a structural barrier and increased dispersion of action potential duration, which may contribute to why patients with structural heart disease are at higher risk for ventricular tachyarrhythmias.
Mesencephalic dopamine neurons ordinarily will not fire faster than about 10/s in response to somatic current injections. However, in response to dendritic excitation, much higher rates are briefly attained. In an analysis of a simplified biophysical model, we suggest a way such high-frequency transient firing may be evoked. Our model represents the neuron as a number of electrically coupled compartments with different natural frequencies, which correspond to the soma and parts of the dendrite. We reduce this model, substituting all the diversity of the compartments that describes real dendritic geometry by a pair of compartments: the slowest, somatic and the fastest, the most distal dendritic one. We have shown that, in the absence of any synaptic stimulation, oscillatory pattern in this model is controlled by the somatic compartment, and, therefore, has a very low frequency. We consider and compare NMDA and AMPA activation applied to the dendritic compartment. Our main result is that activation of the dendritic NMDA receptors evokes oscillations at a much higher frequency. We have also shown that dendritic AMPA activation, by contrast, cannot increase the frequency significantly. The major dynamical question left is how the dendritic frequency can dominate during the application of NMDA. We employ a phenomenon of localization to explain this behavior.
D4 dopamine receptors play an important role in neuropsychiatric disorders involving working memory deficits, such as schizophrenia, attention-deficit/hyperactivity disorder (ADHD) and other disorders of cognition. The impaired methylation has also been documented in different mental diseases. Given a hypothesized role of D4-induced PLM in ion channel functioning [1], we explore in a simple cortical network an influence of this modulation on ability of network to synchronize at gamma frequency and on spike train properties. We model how this mechanism can induce spike trains and modulate its duration and rate, which can mimic short-term memory state. Changes in connectivity within the network influence the effect of DA. Different behavioral characteristics of DRD4 knockout mouse model are believed to be indicative of changes in the striatum (STr), nucleus accumbens (NAc), and prefrontal cortex (PFC), all areas known to be smaller in brains of human ADHD patients. The principal cortical-subcortical neural network for model of schizophrenia has also PFC and NAc/Striatum components. We check how the proposed mechanism can affect transmission in such cortical-subcortical neural network. Results of this study should have significance for further consideration of D4R as novel treatment for ADHD and schizophrenia.
1) Deth, R.C., Kuznetsova A. and Waly, M., 2004, "Attention-related Signaling Activities of the D4 Dopamine Receptor" in Cognitive Neuroscience of Attention, Michael Posner Ed., Guilford Publications, p. 269-282.
I use dynamical systems techniques to model the pheromone response signaling pathway in Saccharomyces cerevisiae yeast. This pathway in yeast is a prototype of regulatory systems that govern response to external stimuli in higher eukaryotic (having a nucleus) organisms. The S. cerevisiae yeast come in two mating types, MATa and MATa. When a yeast cell is exposed to pheromone from yeast of the opposite mating type, a sequence of molecular events happen inside the cell. The components that affect this sequence comprise the pheromone response signaling system.
In my research I focus on about twenty proteins essential to the yeast pheromone response. Each of the proteins involved in the pathway undergoes a state change when the yeast cell is exposed to pheromone. My model is a system of nonlinear coupled differential equations. It describes the production, degradation, phosphorylation (binding of phosphate molecule to a protein), dephosphorylation, and the translocation of proteins as changes in concentrations with respect to time.
The goal is to produce a mathematical model that is simple yet contains enough of the key variables and parameters for accurate predictions of the pathway dynamics. Such a model lends itself to rigorous mathematics (e.g., bifurcation analysis, perturbation theory, stability analysis) and may be used to make predictions that suggest tractable experiments for biologist. There are significant levels of variability in individual cells' ability to respond to pheromone stimulus, thus future work on this project would necessarily include developing a stochastic model.
The yeast Saccharomyces cerevisiae has been used as a model eukaryote for fundamental and applied studies, including stress research. The cellular responses to oxidative stress involve several biological processes, and follow a complex regulation. Understanding which signals trigger the response to the oxidant, which processes are involved, and what is the temporal dynamics of the response is essential. Oxidative stress is related to processes such as ageing, apoptosis, cancer, and the information obtained using yeast cells as a model will certainly help to clarify questions involving higher eukaryotes. In this work we studied the genome-wide response of S. cerevisiae to oxidative stress induced by cumene hydroperoxide (CHP). The changes in gene expression were monitored at the transcriptional level, from a dynamical point of view, spanning a time range of 3 to 120 min after the addition of the oxidant. Affymetrix? Yeast Genome S98 arrays were used for this purpose.
Data obtained in this study was analyzed using several bioinformatics tools and show that oxidative stress dynamics induced by CHP is a complex process. This study also allows an increased resolution about the roles of the genes involved in the oxidative stress response.
Many of the smallest flying insects clap their wings together at the end of each upstroke and fling them apart at the end of each downstroke to augment the lift forces generated during flight. Previous work using rigid wings has shown that at low Reynolds numbers, this mechanism is rather inefficient as large drag forces are produced when the wings are clapped together and pulled apart. In this paper, the immersed boundary was used method to investigate whether or not wing flexibility improves aerodynamic performance during low Reynolds number 'clap and fling.' Our results suggest, for a certain range of flexibilities, that wings with rigid leading edges and flexible trailing edges yield improved aerodynamic performance relative to the rigid wing case. Different wing designs are also shown to yield improved lift generation or improved aerodynamic efficiency.
Graph-theoretic methods are important for the structural analysis of chemical mechanisms. Models that are not capable of replicating experimentally observed behavior can be ruled out by graph-theoretic analysis.
A bipartite graph, used to represent the chemical mechanism connects its structure with the dynamic properties of the corresponding differential equation model. If certain subgraphs are present in the graph then the conventional model can admit instabilities for some values of the system's parameters.
In models with delay or diffusion an instability is usually associated with oscillations or pattern formation respectively. Using the same bipartite graph, the diffusion or delay instabilities can be characterized also in terms of the structure of the chemical mechanism.
Molecular motors, such as kinesin, myosin, or dynein, convert chemical energy into mechanical energy, by hydrolyzing ATP. The mechanical energy is used for moving in discrete steps on the cytoskeleton, and carrying a molecular load. With sufficient resolution, single molecule recordings of motor steps appear as a stochastic sequence of dwells, resembling a staircase. Here, we developed maximum likelihood algorithms that separate the dwell time sequence from noise, estimate the mean and variance of the step size between consecutive substrate positions, and estimate the rate constants of conformation and step transitions of the molecular motor. We model the motor with an infinite but periodic Markov model, reduced to a small model reflecting the periodic chemistry of each step. The dwell sequence is extracted using maximum likelihood dynamic programming algorithms. The kinetics are estimated from the dwell time sequence by numerical maximization of the likelihood fun! ction for discrete time Markov models. The algorithm can fit models with arbitrary chemistry and allows global fitting across stationary and nonstationary experimental conditions, and user-defined constraints on rate constants. Our results show that when the measurement noise is within the nanometer range, steps as small as 8 nm can be analyzed. The algorithm is implemented in the free QuB software (www.qub.buffalo.edu).
In order to study how the convergence of many variable neurons on a single target can sharpen timing information, we investigate the limit as the number of input neurons and the number of incoming spikes required to fire the target both get large with the ratio fixed. We prove that the standard deviation of the firing time of the target cell goes to zero in this limit and we derive the asymptotic forms of the density and the standard deviation near the limit. We use the theorems to understand the behavior of octopus cells in the mammalian cochlear nucleus.
In this paper, classical epidemiology models are modified to incorporate the effects of control stages. Effects of isolation and qua-ranting the traced contacts to wipe out the disease from the population are highlighted. Three models with different contact tracing functions for new emerging diseases are described and studied. A model with general contact tracing rate is also framed. Threshold, disease free equilibrium and its global stability criteria for the models are explored.
Analysis of parameters representing quarantine and isolation efforts, on the threshold is done. Data from the SARS outbreak in Hong Kong is used in numerical simulations to illustrate the results. Sensitivity and uncertainty analysis on the traced reproduction number is performed. Cost analysis is detailed and carried out for these control measures. We determine the effect of different interventions on the traced reproduction number and estimate requirements to achieve elimination of the infectious disease. We find that the strategy of tracing and quarantining contacts of diagnosed cases can be very successful in reducing transmission, but large scale contact tracing efforts might prove economically prohibitive.
Morphogen gradients of Wingless are involved in patterning the Drosophila embryo and imaginal discs. The dpp and wg genes are expressed in the dorsal and ventral regions of the drosophila leg imaginal disc. A key factor in maintaining the non-overlapping regions of gene expression is the inhibition of dpp transcription by Wg signaling. The recent experimental studies of ventral leg disc cell showed the Arm/dTcf complexes activate the wg expression. In addition, the Arm/dTcf complexes cooperate with Brinker (Brk) to inhibit dpp transcription. Furthermore, ectopic expression of Arm or dTcf causes loss of dpp repression. To validate the experimental results and to determine the mechanism that leads to dpp repression, we developed an intra-cellular mathematical model to investigate the protein-protein and protein-DNA interactions. Our model showed that dpp repression occurred in response to Wg singaling only when Brk interacts with Arm to form a repressing complex R! 1TABR2 in the presence of non-repressing complexes. In the wing disc, Wingless signaling is affected by binding of the ligand to the glycosaminoglycan (GAG) chains of proteoglycans. GPI anchored heparan sulfate proteoglycans (HSPG) such as glypicans, antagonize wingless signaling by trapping it in the extracellular matrix. Syndecan (Sdc) is another kind of proteoglycan that is decorated with sugar chains. Current study shows that Drosophila Syndecan also modulates Wingless signaling but by a different mechanism. The one syndecan gene is transcribed to multiple isoforms with specific tissue distribution. Experimental results suggest that Syndecan antagonizes Wingless signaling by promoting the internalization of the secreted Wingless protein. Our extra-cellular mathematical model incorporates diffusion, reaction, internalization of ligand, ligand-receptor and ligand-nonreceptor to investigate such system. We found that when Sdc is overexpressed, loss of Wg signal occur! s and the amount of internal ligand complex increases.
We develop a deterministic mathematical model to describe reactivation of latent virus by chemical inducers. This model is applied to the reactivation of latent KSHV in BCBL-1 cell cultures with butyrate as the inducing agent. Parameters for the model are first estimated from known properties of the exponentially growing, uninduced cell cultures. The model is then extended to describe chemically induced KSHV reactivation in latently infected BCBL-1 cells. Additional parameters that are necessary to describe induction are determined from fits to experimental data from the literature. Our model provides good agreement with two independent sets of experimental data.
To examine the network-level encoding used by hippocampus to achieve its real-time representations of episodic information, we have used a 96-channel array to simultaneously record the activity patterns of as many as 260 individual neurons in the mouse hippocampus during various startling episodes. We find that the mnemonic startling episodes triggered firing changes in a set of CA1 neurons in both startle-type and environment-dependent manners. Pattern classification methods reveal that these firing changes form distinct ensemble representations in a low-dimensional encoding subspace. Application of a sliding window technique further enabled us to reliably capture not only the temporal dynamics of real-time network encoding but also postevent processing of newly formed ensemble traces. Our analyses revealed that the network-encoding power is derived from a set of functional coding units, termed neural cliques, in the CA1 network. The individual neurons within neural cliques exhibit ''collective cospiking'' dynamics that allow the neural clique to overcome the response variability of its members and to achieve real-time encoding robustness. Conversion of activation patterns of these coding unit assemblies into a set of real-time digital codes permits concise representation and categorization of discrete behavioral episodes.
We study the problem of closure of a hierarchy of coupled integro-differential equations describing the dynamics of the product densities m[k],k=1,2,..., arising from the deterministic approximation of a locally regulated, spatial birth-and-death point process, with rules motivated by plant population dynamics. In our approach, we truncate at order k=2 and close the product density m[3] through maximisation of the entropy functional of the process, restricted to a small (but not infinitesimal) observation window A, in such a way that the likelihood of observing configurations involving more than three points in A is negligible and A captures the scale of spatial correlations of third order. The maximization is carried out subject to the constraints of normalization and fixed product densities m[1] and m[2], that are given by the hierarchy. We thereby obtain an implicit, Kirkwood-type closure for m[3] , coupled to an ancillary integral equation for the domain of triplet correlations A. This new closure is enhanced with previously unnoticed correction terms, which vanish for the homogenous Poisson process. They are shown to become important as the process deviates from complete spatial randomness. Numerical solutions of the product density dynamics, based on this new closure, are in remarkably good agreement with the equilibrium values from individual-based simulations, both for aggregated and segregated spatial patterns.
Angiogenesis, the formation of new blood vessels from pre-existing vessels, is a normal process in growth and development, as well as in wound healing. Recent experimental work showed the importance of the capillary sprouts in the formation of a new micro-vasculature in a wound. Hence, we focus our attention on sprouting angiogenesis and blood microcirculation in a healing wound.
Our interest is to model the plasma flow into a cylindrical axisymmetric capillary sprout with permeable walls and non-uniform diameter. According with the form of the capillary, we consider three main cases: first, the tube is narrow enough to exclude red blood cells, second, the red blood cells are stacked inside the sprout forming a kind of porous media, and third, where the capillary is wide enough for red blood cells to be squashed by the pressure-driven plasma flow, forming a poroelastic medium.
We assume that the fluid pressure is greater than the pressure outside the channel, accounting also for osmotic pressure differences. Thus, we expect to have only outward flow through the walls and also, we define the outward flux per unit area using Starling's law of membrane filtration.
The plasma flux into the sprout as a function of the wall's permeability and sprout shape is determined. In particular, the flux through the sprout increases with sprout length, because of the increased area over which plasma can leak through the wall. The results are compared with in vivo experimental observations.
The plasma inflow is determined by the microcirculation of blood through the main network. Our interest is to compute the pressures and fluxes of a given capillary network and to study its stability. We consider the most simple network which include all the main aspects of the microcirculation, network introduced by Carr[1]. The rheological properties of the blood are modeled using the Pries-Secomb[2,3] in vivo determined relations. The stability of the system, for various parameters, is discussed.
References:
- R. T. Carr, J. B. Geddes and F. Wu, "Oscillations in a simple microvascular" network, Annals of Biomedical Engineering, Vol. 33, 6, 764-771(8), 2005.
- A. R. Pries, T. W. Secomb, P. Gaehtgens and J. F. Gross, "Blood flow in microvascular networks. Experiments and simulation", Circulation Research, Vol. 67, 826-834, 1990.
- A. R. Pries, T. W. Secomb, T. Gessner, M. B. Sperandio, J. F. Gross and P. Gaehtgens, "Resistance to blood flow in microvessels in vivo", Circulation Research, Vol. 75, 904-915, 1994.
We introduce an integrodifference equation model to study the spatial spread of epidemics through populations with overlapping and non-overlapping epidemiological generations. Our focus is on the existence of travelling wave solutions and their minimum asymptotic speed of propagation c*. We contrast the results here with similar work carried out in the context of ecological invasions. We illustrate the theoretical results numerically in the context of SI (susceptible-infected) and SIS (susceptible-infected-susceptible) epidemic models.
Even though much is known about the biophysics, anatomy and physiology of basal ganglia networks in Parkinson's disease, there is no clear understanding of how all these properties are connected to each other and what is the origin of the motor symptoms, including parkinsonian tremor. We present the development of biophysically based mathematical model of the basal ganglia-thalamocortical circuits of the brain, which are involved in the generation of pathological tremor. The model suggests the mechanism of how the cellular properties of pallidal and subthalamic neurons may give rise to tremor oscillations in the presence of thalamocortical feedback, thus providing insights into pathophysiology of parkinsonian tremor.
An averaged ODE model is presented to explore the effects of macrophage infiltration in tumours. The tumour is portrayed as an ensemble of three interacting populations: proliferating cells, quiescent cells and necrotic material. The process of vascularization is also taken into account and macrophages are chemotactically attracted into the tumour. Different scenarios are considered: either the quiescent cells or the necrotic material could promote the chemotactic recruitment and macrophages can be engineered to have therapeutic effects such as cytotoxicity or inhibition of angiogenesis. A bifurcation analysis is presented exposing qualitative differences in the dynamics of the possible submodels, e.g. the onset of oscillatory behaviour and the appearence of multiple stable stationary solutions. Special attention is payed to those contrasting behaviours which seem more relevant when devising future laboratory experiments.
FtsZ, a bacterial homolog of tubulin, forms single-stranded polymers that assemble cooperatively. Models for cooperative polymerization traditionally require polymers to be multistranded, suggesting that new models are now needed. Potentially, cooperativity might emerge if a subunit changes from a low to a high affinity conformation when in contact with adjacent subunits in the polymer. Computer programs that model chemical reactions could determine whether such models for single-stranded polymers produce the lags and critical concentrations indicative of cooperativity. However, conventional programs are difficult to apply to polymers because an unlimited number of different polymer species may exist, whereas these programs require discretely defined species. Additionally, conventional programs cannot efficiently track polymer lengths, the spatial relationships between subunits within a filament or make rate constants dependent upon them. We have developed a computer program that can track the behavior of 10,000 subunits based on probabilities of reaction per unit time. The modeling strategy uses a stochastic reaction event generator. To validate the algorithm, we have successfully simulated mathematically solvable polymerization schemes, including isodesmic polymerization and the model for double-stranded polymerization used to fit FtsZ data1. Our program is robust, eliminating feedback loops that can occur in deterministic models when all long polymers are combined into one chemical class. Additionally, our program can track changing polymer length distributions over time. We are currently extending the program to determine whether cooperativity in single-stranded FtsZ polymers can be explained by models of polymerization that include nucleotide hydrolysis, or in which a subunit's reaction rates are tied to the chemical states of its neighbors.
Joint work with Dr. Laura Romberg, Oberlin College Dept. Biology.
Dendrites of many nerve cells are complex, branching structures that receive and process thousands of synaptic inputs from other neurons. Loci for receiving inputs are served by dendritic spines that are equipped with excitable channels. Here we develop a mathematical model of branched dendritic tree based upon a generalisation of the analytically tractable Spike-Diffuse-Spike model. The active membrane dynamics of spines are modelled by an integrate-and-fire process. The spines are assumed to be discretely distributed along a passive branched dendritic structure. The model supports saltatory travelling wave propagation and wave scattering amongst the dendritic branches. This model is ideally suited for the study of spatio-temporal filtering properties and neural response to different patterns of synaptic input.
An important problem in computational biology is the modeling of several types of networks, ranging from gene regulatory networks and metabolic networks to neural response networks. In [LS], Laubenbacher and Stigler presented an algorithm that takes as input time series of system measurements, including certain perturbation time series, and provides as output a discrete dynamical system over a finite field. Since functions over finite fields can always be represented by polynomial functions, one can use tools from computational algebra for this purpose. The key step in the algorithm is an interpolation step, which leads to a model that fits the given data set exactly. Due to the fact that biological data sets tend to contain noise, the algorithm leads to over-fitting.
Here we present a genetic algorithm that optimizes the model produced by the Laubenbacher-Stigler algorithm between model complexity and data fit. This algorithm too uses tools from computational algebra in order to provide a computationally simple description of the mutation rules.
We describe an application of the combined algorithm in a computational neuroscience project, that in collaboration with a neuroscience research group in the Rutgers Psychology Department, has as ultimate goal to apply our modeling techniques to fMRI data collected from spinal cord patients in order to study possible venues for pain management.
[LS] Laubenbacher, R. and B. Stigler, A computational algebra approach to the reverse-engineering of gene regulatory networks, J. Theor. Biol. 229 (2004) 523-537.
Co-authors: Lan Ma(1), John Jeremy Rice(1), Wenwei Hu(2), Arnold Levine(2) and Gustavo Stolovitzky(1)
(1) IBM T.J. Watson Research Center, Yorktown Heights, NY
(2) The Cancer Institute of New Jersey, Robert Wood Johnson School of Medicine, New Brunswick, NJ
The tumor suppressor p53 is critical to ensure genomic stability. In single cells, the oscillatory p53 response to ionizing radiation (IR), which induces double stranded breaks (DSBs), is "digital," in that the number of oscillations rather than the amplitude shows dependence on the radiation dose. We present a model of single cell p53 dynamics in response to ionizing radiation. In our model, DSB sites interact with a limited pool of DNA repair proteins, forming DSB-protein complexes at DNA damage foci. Both the initial number of DSBs and the DNA repair process are modeled stochastically. The model assumes that the persisting complexes are sensed by ataxia telangiectasia mutated (ATM) kinase, which transduces in an ON/OFF manner the DNA damage signal to the downstream negative feedback oscillator consisting of p53 and its negative regulator Mdm2, a transcriptional target of p53. Our model exhibits coordinated oscillations of p53 and Mdm2 upon IR stimulation, with a stochastic number of oscillations whose mean increases with IR dosage, in agreement with the observed response of p53 to DNA-damage in single-cell experiments.
The study of the motility of cilia and flagella is of great importance in biology and medicine. Mathematical modeling to simulate the internal mechanism of cilia and flagella has attracted many researchers. But no satisfactory work has been done for multiciliary beating. We present an integrative computational model of multiciliary beating. The axoneme is modeled by Dillon-Fauci approach which was first used in [1] in 2000. Dillon-Fauci axoneme model is simple but very successful in modeling the cilia and sperm motility (see [1,2]). This model, based upon the immersed boundary method (Peskin), couples the internal force generation of the molecular motors through the passive elastic structure with the external fluid mechanics governed by the Navier-Stokes equations. In our model, the multiciliary configuration by the immersed boundary method does not cause much extra cost in computation. In our numerical results, we will show how a single cilium interacts with its n! eighboring cilia, how viscosity effects its beating frequency, how metachronal wave, synchronization phenomena are genereated. The challenging problems wherein will also be discussed. At last we present the computer simulations for chlamydomonas swimming and muco-ciliary interactions as applications of our model. This project is a joint work with Professor Lisa Fauci, Tulane University and Robert Dillon, Washington State University.
References:
- Robert H. Dillon and Lisa J. Fauci, "An Integrative Model of Internal Axoneme Mechanics and External Fluid Dynamics in Ciliary Beating", J. Theor. Biol. v. 207, p. 415--430, 2000.
- Robert H. Dillon, Lisa J. Fauci, and Charlotte Omoto, Mathematical modeling of axoneme mechanics and fluid dynamics in ciliary and sperm motility. Dyn. Contin. Discrete Impuls. Syst. Ser. A, Math. Anal., 10(5):745-757, 2003. Progress in partial differential equations (Pullman, WA, 2002).
