2008 Workshop for Young Researchers in Mathematical Biology (WYRMB): Abstracts and Lecture Materials
Plenary Speakers
Vaginal microbicides, designed to prevent HIV infection in women, are one of the most promising biomedical interventions. Clinical trials of second-generation microbicides have begun; if shown to be effective they could be licensed within 5-10 years. Since these microbicides contain antiretrovirals (ARVs) they could be highly effective. However there is concern that, if used by HIV-positive women, ARV resistance may evolve. By analyzing a novel mathematical model we find that adherence could have both beneficial and detrimental effects on trial outcomes. Most importantly we show that planned trial designs could mask resistance risks and therefore enable high-risk microbicides to pass clinical testing. We then parameterize a transmission model using epidemiological, clinical and behavioral data to predict the consequences of wide-scale usage of high-risk microbicides in a heterosexual population. Surprisingly, we show that reducing a participant's risk of resistance during a trial could lead to unexpectedly high rates of resistance afterwards when microbicides are used in public health interventions. We also find that, paradoxically, although microbicides will be used by women to protect themselves against infection they could provide greater benefit to men. More infections in men, than women, will be prevented if there is a high probability that ARVs are systemically absorbed, microbicides are less than ~50% effective, and/or adherence is less than ~60%. Men will always benefit more than women in terms of infections prevented per resistant case; but this advantage decreases as the relative fitness of drug-resistant strains increases. Interventions that use ARV-based microbicides could have surprising consequences.
Understanding the dynamics of human body weight change has important consequences for conditions such as obesity, starvation, and wasting syndromes. Changes of body weight are known to result from imbalances between the energy derived from food and the energy expended to maintain life and perform physical work. However, quantifying this relationship has proved difficult, in part because the body is composed of multiple components and weight change results from alterations of body composition (i.e., fat versus lean mass). Here, I will describe recent modeling efforts that provide a general description of how body weight will change over time by tracking the flux balances of the macronutrients fat, protein, and carbohydrates. The model provides an explanation of the recent obesity epidemic, dispels some myths about weight loss and weight gain, and can explain the phenomenon of "catch-up" growth in children.
Among all the regulators of cellular signalling, Calcium ranks as one of the most important. In many cell types Calcium signalling is associated with the generation of waves. In this talk I will present a bidomain fire-diffuse-fire model that facilitates mathematical analysis of such propagating waves in living cells. Modelling Calcium release as a threshold process allows the explicit construction of travelling wave solutions to probe the dependence of Calcium wave speed on physiologically important parameters. A linear stability analysis identifies the newly discovered 'Tango waves' as a dynamic instability of propagating pulse solutions. Furthermore, I will use the bidomain model to shed some light on recently discovered sensitisation waves in ventricular myocytes. Moving next to atrial myocytes, I will describe the development of a computationally efficient framework to investigate stochastic Calcium release and propagation in a realistic three-dimensional cellular geometry.
In a regression or classification problem, one often has many potential predictors (independent variables), and these predictors may interact with each other to exert non-additive effects. I will present a Bayesian approach to search for these interactions. We were motivated by the epistasis detection problem in population-based genetic association studies, i.e., to detect interactions among multiple genetic defects (mutations) that may be causal to a specific complex disease. Aided with MCMC sampling techniques, our Bayesian method can efficiently detect interactions among many thousands of genetic markers.
A related problem is to find patterns or associations in text documents, such as the "market basket" problem, in which one attempts to mine association rules among the items in a supermarket based on customers' transaction records. For this problem, we formulate our goal as to find subsets of words that tend to co-occur in a sentence. We call the set of co-occurring words (not necessarily orderly) a "theme" or a "module". I will illustrate a simple generative model and the EM algorithm for inferring the themes. I will demonstrate its applications in a few examples including an analysis of chinese medicine prescriptions and an analysis of a chinese novel.
Vascular development and homeostasis are underpinned by two fundamental features: the generation of new vessels to meet the metabolic demands of under-perfused regions and the elimination of vessels that do not sustain flow. I will present a multiscale model of vascular tissue growth that combines blood flow, angiogenesis, vascular remodelling and the subcellular and tissue scale dynamics of multiple cell populations. Simulations show that vessel pruning is highly sensitive to the pressure drop across a vascular network, the degree of pruning increasing as the pressure drop increases. In the model, low tissue oxygen levels alter the internal dynamics of normal cells, causing them to release Vascular Endothelial Growth Factor (VEGF), which stimulates angiogenic sprouting. Consequently, the level of blood oxygenation regulates the extent of angiogenesis, with higher oxygenation leading to fewer vessels. Simulations show that network remodelling (and de novo network formation) is best achieved via an appropriate balance between pruning and angiogenesis. An important factor is the strength of endothelial tip cell chemotaxis in response to VEGF. When a cluster of tumour cells is introduced into normal tissue, as the tumour grows hypoxic regions form, producing high levels of VEGF that stimulate angiogenesis and cause the vascular density to exceed that for normal tissue. If the original vessel network is sufficiently sparse then the tumour may remain localised near its parent vessel until sufficient new vessels bridge the gap to an adjacent vessel. This can lead to metastable periods, during which the tumour burden is approximately constant, followed by periods of rapid growth.
The hippocampal formation is crucial for remembering episodes in one's life, and evidence suggests that synchronous activity throughout the hippocampus is essential for the mnemonic functions of this brain structure. We have studied the mechanisms of synchronization using electrophysiological and computational methods. More recently, we have exploited methods for introducing real-time control in cellular electrophysiology. These techniques allow us to "knock in" virtual ion channels that can be controlled with great mathematical precision, and to immerse biological neurons in real-time, virtual neuronal networks. These manipulations allow us to test computationally-based hypotheses in living cells. From this work, I will discuss which properties of single cells and neural networks seem crucial for coherent activity in the hippocampal formation. If time permits, I will also discuss work on the consequences of precise spike timing in neuronal network function.
Selected related papers:
- Netoff TI, Banks MI, Dorval AD, Acker CD, Haas JS, Kopell N, and White JA (2005) Synchronization in hybrid neuronal networks of the hippocampal formation. Journal of Neurophysiology 93: 1197-1208.
- Dorval AD and White JA (2005) Channel noise is essential for perithreshold oscillations in entorhinal stellate neurons. Journal of Neuroscience 25: 10025-10028.
- Zhou YD, Acker CD, Netoff TI, Sen K, and White JA (2005) Increasing calcium transients by broadening postsynaptic action potentials enhances timing-dependent synaptic depression. Proceedings of the National Academy of Sciences USA 102: 19121-19125.
- Keck T, Lillis KP, Zhou YD, and White JA (2008) Frequency-dependent glycinergic inhibition modulates plasticity in hippocampus. Journal of Neuroscience 28: 7359-7369.
- Fernandez FR and White JA (2008) Artificial synaptic conductances reduce subthreshold oscillations and periodic firing in stellate cells of the entorhinal cortex. Journal of Neuroscience 28: 3790-3803.
MBI Postdoc Presentations
It is known that adipose tissue, formed by fat cells, is a regulatory organ for energy homeostasis and responsible for thermogenesis. Recent discoveries show that the tissue is also an endocrine organ. Pathology of the adipose tissue is interrelated with obesity, diabetes, insulin resistance, and atherosclerosis. Understanding the mechanism of the tissue, especially at the transcriptional level, is thus critical for developing clinical methods for the treatment. In this article a mathematical model is developed and analyzed. The analysis suggests that the decision is made among different states of cell fate, such as differentiation, proliferation and quiescence, depends mainly on the inhibition of NF-$\kappa$B. The model provides a quantitative measure for prediction of the fate of a preadipocyte based on some parameter values. This result along with experimental justification are the main points addressed in the article.
Controlling inflammation has become a key focal point in the treatment of critically ill patients and the identification of proper biological targets and their specific manipulations is necessary. One means of achieving this end involves formulating a strategy for delivering therapies, in the correct amount, at the right time. A tool that can help determine this complex dose regimen is Nonlinear Model Predictive Control (NMPC). We apply NMPC in the setting of a 4-dimensional ordinary differential equations model of the acute inflammatory response to pathogen due to Reynolds and colleagues (JTB, 2006, v. 242).
A patient population of 1000 randomly generated patients was considered. Of these virtual patients, 620 were recognized as being ill enough to receive treatment, as defined by levels of phagocytic cells exceeding a predetermined threshold. The algorithm identified individualized therapeutic strategies that enabled approximately 20-40% more patients to be successfully restored to homeostasis (damage reduced to baseline) as compared to other treatment strategies (e.g. non-individualized therapies or placebo). Current and future work will also be discussed.
A fundamental question in invasion biology and wildlife management concerns the most effective methods to limit the spread of and to control a harmful population. This question has been analyzed in a variety of theoretical frameworks.
The hypothesis we investigate is whether optimal control theory applied to an aggregated, analytic model can be used to effectively control a harmful species modeled by an IBM, or whether interactions between individuals, their spatial distribution, and landscape heterogeneities limit the effectiveness of the control methods derived from the analytic model. If optimal policies derived from an aggregated model are determined to be generally effective in managing a population modeled with considerably greater complexity, this provides evidence that optimal management strategies may be relatively insensitive to the details of individual behavior and the associated impacts on population response. Even if this does not hold in general, any cases for which it does hold can provide potentially useful "rules-of-thumb" for population management.
Motivated by the aim of modelling the mechanical behaviour of biological gels (such as collagen gels) which have a fibrous microstructure, we consider the flow of a two-dimensional film of incompressible, transversely isotropic viscous fluid. Two regimes are considered: extensional flow, and squeezing flow of the fluid between two rigid plates. Neglecting inertia, gravity and surface tension, in each regime leading-order equations are derived from the full flow problem using a thin-film approximation. Special cases in which the solution may be determined explicitly are considered and the physical interpretation of the results is discussed.
Glioma cells, already at an early growth of the tumor, tend to migrate from the primary tumor into the surrounding tissue in different patterns. For example, cells from U87 cell-line have been shown to disperse in a radially symmetric fashion from a spherical tumor, whereas mutant U87DeltaEGFR cells form branching patterns similar spokes from a hub. Several models have been suggested to explain the different modes of migration, but none of them, so far, has explored the important role of cell-cell adhesion. The present paper develops a mathematical model which includes the role of adhesion and provides an explanation for the various patterns of cell migration. It is shown that, depending on critical adhesion and chemotactic parameters, the migration patterns exhibit a gradual shift from branching to dispersion, as observed in experiments.
Co-authors: Avner Friedman, Sean Lawler, Michal O. Nowicki and E. Antonio Chiocca
In the study of calcium waves, the endoplasmic rheticulum (ER) and its IP3 mediated calcium induced calcium release (CICR) mechanism have dominated the spotlight. However, more than the ER can display calcium dependent excitability. The oft-times overlooked (in regards to Ca2+ signaling) mitochondria also exhibit CICR via the permeability transition pore (PTP). In our work, we extend the well-known Fall-Keizer mathematical model for the Ca2+ dynamics between the ER, the cytosol, and the mitochrondria to include dynamics due to the PTP. Three key components come into play: non-steadystate mitochrondrial proton concentration, the permeability transition pore, and a weak-acid term that is essential for robust mitochondrial homeostasis. We show that, with these simple additions, the model exhibits mitochondrial CICR (mCICR). Furthermore, we extend the model from a single point to a spatially extended system and find traveling waves of both \Ca and potential, as found by Ichas et al (1997). Additionally, the model also captures the pore "popping" event whereby the pore opens and remains in a high-conductance state, which ultimately leads to apoptosis or necrosis due to calcium cytotoxicity.
I will present a biologically-based mathematical model that accounts for several features of the human sleep/wake cycle. These features include the timing of sleep and wakefulness under normal and sleep-deprived conditions, ultradian rhythms, and the circadian dependence of several sleep variables such as the total sleep time and sleep-onset rapid eye movement sleep (SOREMS). Additionally, if the input from the neurotransmitter orexin is removed, the system exhibits more frequent switching between sleep and wakefulness, consistent with the sleep disorder narcolepsy. The model demonstrates how each of these features depend on interactions between a circadian pacemaker and a sleep homeostat and provides a biological basis for the two-process model for sleep regulation(Borbely, 1982; Daan et al., 1984). The model is based on previous "flip-flop" models for sleep/wake (Chou et al. 2001) and REM/NREM (Lu et al. 2006) and we explore whether the neuronal components of these flip-flop models, with the addition of a sleep-homeostatic process, are sufficient to account for the features of the sleep/wake cycle listed above. The model is minimal in the sense that, besides the sleep homeostat and constant cortical drives, the model includes only those nuclei described in the flip/flop models. However, in order to account for certain features of the ultradian rhythms, we found it necessary to add an additional hypothesis about the connections.
Wound healing represents a well-orchestrated reparative response that is induced by injuries. Angiogenesis, the formation of blood vessels from existing vasculature, plays a central role in wound healing. In this talk, I will present a mathematical model that addresses the role of tissue oxygen tension in cutaneous wound healing. Key components of the developed model include capillary tips, capillary sprouts, fibroblasts, inflammatory cells, chemoattractants, oxygen, and the extracellular matrix. The model consists of a system of nonlinear partial differential equations describing the interactions in space and time of the above variables. The simulated results agree with the reported literature on the biology of wound healing. The proposed model represents a useful tool to analyze strategies for improved healing and can be used to generate novel hypotheses for experimental testing.
DNA methylation has been shown to play an important role in the silencing of tumor suppressor genes in various tumor types. As it is desirable to have a system-wide understanding of the methylation changes that occur in tumors, a differential methylation hybridization (DMH) protocol has been developed (Yan, et al. 2008) and it can simultaneously assay the methylation status of all known CpG islands. As is common with all microarray technologies, a large percentage of the signal obtained from the array can be attributed to various measurable and unmeasurable confounding factors unrelated to the biological question at hand. In order to correct the bias due to the noise, we have implemented a linear regression model for assessing significance of signal. In particular, we focus on finding the CpG islands that are differentially methylated understand two conditions, and the highly methylated CpG islands that can identify more aggressive tumors. We have been able to correctly identify some known methylated and unmethylated genomic regions.
GM-CSF is a drug that enhances the ability of macrophages to present antigen and initiate immune response. GM-CSF also stimulates monocytes to secrete soluble VEGF receptor (sVEGFR-1) which binds to and inactivates VEGF. Eubank and colleagues in a recent work discovered that GM-CSF treatment locally in murine breast cancers reduced tumor growth and metastasis; moreover, GM-CSF lowered oxygen level and reduced blood vessel density within the tumor. We developed a mathematical model which addresses the effect of GM-CSF on the growth of breast cancer in mice. The model takes into account the experimentally established interactions among cancer cells, macrophages, endothelial cells, VEGF and M-CSF. The model simulates the growth of tumor as a function of the local GM-CSF dose. The model simulations were validated against in vivo data and show a good fit with experimental results. We used the model to compare the efficacy of different dosing protocols of injection of GM-CSF, as well as to suggest new hypotheses for slowing the progression of breast tumor. For example, the model suggests that injecting the drug daily, twice or 3 times a week are comparably effective. In contrast, reducing or over increasing the frequency of dosing is counterproductive. We hope to further refine the model in future work, by including the interactions of macrophages and other immune cells, fibroblasts and cytokines that communicate between the tumor and its micro-environment.
Co-authors: Tim D. Eubank, Ryan D. Roberts, Clay B. Marsh, and Avner Friedman
Poster Presentations
Epidemiological models can be used to make predictions of the future time course of an outbreak. In order to have confidence in such predictions, we must have reliable estimates of model parameters. Unfortunately, as with most biological data, epidemic data pose a number of challenges that need to be overcome before they can be relied on to provide parameter estimates. We take an in-depth look at parameter estimation for the Susceptible-Infective-Recovered (SIR) model when presented with different scenarios of data availability. We use asymptotic statistical theory, together with sensitivity analysis, to quantify uncertainty of the estimated parameters. We examine the importance of specific data points at certain times throughout the epidemic to find when to collect data to best obtain information for parameter estimation. We conclude that there are optimal time intervals during an outbreak when increased data collection results in better estimations. Epidemiological data is also subject to under-reporting; oftentimes, not all cases of infection during an epidemic are reported to data collectors. We investigate how under-reporting impacts parameter estimation and show that failure to account for under-reporting yields unrealistically low values of the basic reproductive number resulting in overly optimistic predictions of epidemics.
The main goal of this work is to determine an optimal treatment strategy using interferon and ribavirin, through mathematical modeling. We formulate a mathematical model using a system of ordinary differential equations, which describes the interaction of target cells, infected cells, infectious virions, non-infectious virions and the two drugs. We analyze the model for steady states and do a stability analysis. We solve an optimal control problem with an objective functional that maximizes the target cell population. Finally we derive optimal treatment strategy and then solve it numerically.
A mathematical model to study the effects of HER2 overexpression on cell proliferation in breast cancer is presented. The model illustrates the proliferative behavior of cells as a function of HER2 and EGFR receptors numbers. This mathematical model comprises kinetic equations describing the cell surface binding of EGF growth factor to EGFR and HER2 receptors, coupled to a model for the dependence of cell proliferation rate on growth factor receptors binding.
One of the first lines of defense against tuberculosis is the bacille Calmette-Guérin (BCG) vaccine. Although BCG is an old vaccine, studies of its protective efficacy have had widely divergent results. In addition, the vaccine can interfere with the diagnosis of latent TB infection. Due to these qualities, the international community has adopted varying policies on the vaccine's use.
We introduce a dynamical system model for tuberculosis dynamics at the population level in order to examine the conditions which justify discontinuing BCG vaccination. Whenever possible, model parameters are defined in terms of control indicators estimated by the World Health Organization. The model is analyzed via numerical experiments, sensitivity analysis and dynamical systems theory. The model is fit to TB data for ten countries with varying TB burdens. In addition to the role of vaccine, we can consider alternative intervention strategies on a country-specific basis.
The two-stage clonal expansion (TSCE) model of carcinogenesis is a mechanistic model that has been used to successfully explain epidemiologic patterns of lung cancer risk. We derive a modified two-stage clonal expansion model to explain the excess lung cancer risk patterns in data from the Mayo Lung Project (MLP), a randomized clinical trial designed to assess the efficacy of chest Xray screening in reducing lung cancer mortality. More traditional risk estimation methods have low power to detect levels of excess cancer risk such as observed in the MLP. An advantage of our model-based framework is that it allows us to incorporate correlates of excess lung cancer risk such as age and screen frequency into a likelihood-based estimation of the model parameters. This methodology greatly increases the power to detect a difference in excess lung cancer risk attributable to screening when one exists, and when the relationship between age, screen frequency, and excess lung cancer risk are consistent with the model-based predictions. In this study, we describe our mathematical formulation of excess lung cancer risk based on the TSCE model for a population of high-risk smokers. We estimate the parameters of our model directly from long-term follow-up data in the Mayo Lung Project, using a Bayesian data augmentation procedure combined with maximum likelihood estimation, and examine the significance of our parameter estimates. Finally, we examine the goodness of fit of our simulation model of lung cancer onset and tumor progression both with and without the excess risk parameters to the long-term follow-up data from the MLP.
Work done in collaboration with Marek Kimmel.
In early Drosophila development, one of the first specifications made in determining cell location is the formation of the anterior-posterior axis, and this event is carried out with incredible precision, despite many variations of the genetic regulatory network and environmental fluctuations. My focus is on the first step in this process, where the protein, Hunchback (Hb), forms a very sharp boundary at the midpoint of the embryo, despite receiving noisy information from its known upstream regulator, Bicoid (Bcd). Using systems of partial differential equations to model the gene network, we seek an answer to the following inverse problem: given the robust Hb boundary, what conditions does this place on the regulatory network? While a simple answer does not likely exist, we provide examples of such networks, as the result of both numeric searches and analytic techniques, and an interpretation of the conditions in terms of the biological system.
Since its introduction into eastern North America in the 1940s, the eastern population of house finches (Carpodacus mexicanus) has become partially migratory, unlike its on-migratory source population in southern California (Able and Belthoff, 1998; Belthoff and Gauthreaux, 1991). The infectious disease mycoplasmal conjunctivitis (pathogen Mycoplasma gallisepticum or MG), which has been monitored in the house finch population since its appearance around 1993 (Dhondt et al., 1998), may induce higher mortality rates among populations in more northerly latitudes relative to more southerly populations. Here we investigate the potential impact of this differential disease mortality on the migratory structure of the eastern house finch population using an epidemic modeling approach (forced ODEs). Analytical and computational results suggest the ongoing MG epidemic in the eastern house finch could lead to increases in the percentage of and the total number of migrating individuals in a population despite overall population declines, assuming relatively high winter mortality rates in the north eastern part of their range. These results also suggest that empirical evidence of such a change in migratory structure could be most noticeable in northerly inland populations that showed significant declines following the initial outbreak of MG in the east.
Synchronization of coupled oscillators has been an area of interest for a long time. Depending on the intrinsic properties of the oscillators and the coupling between them, it is possible to observe different activity patterns in the network. Here, we use the phase resetting theory to uncover the synchronization properties of coupled neural oscillators. Existence and stability criteria for harmonic locking modes were derived for two reciprocally pulse coupled oscillators based on their first and second order phase resetting curves (PRCs). These methods were then tested using two reciprocally inhibitory Wang-Buzsaki model neurons.PRCs are generated in an open loop configuration with the assumption that the synaptic inputs received in the closed loop circuit remain similar to those used to generate the PRCs. We assume a firing pattern and use it to produce a map, which can then be linearized for a stability analysis. A periodicity criteria is derived in terms of the stimulus and response intervals for each neuron for every cycle which results in the iterative map. A stability analysis of the linearized map provides a single eigenvalue for the map, provided the second order resetting of all but the last input in a cycle is disregarded. Previously, Ermentrout derived existence and stability criteria for n:m locking assuming weak coupling using averaging theory. The methods presented here do not require the coupling to be weak and are easier to implement and apply to real neurons since only the PRCs are required. Both methods agree in the weak coupling regime, but the new method shows good agreement with the observed modes from the simulated network even in strong coupling regimes. Our analysis show that the information obtained from PRCs is sufficient to determine synchronization properties of the neural network.
We have investigated how noise propagation in biological reaction networks affects system sensitivities. If the numbers of molecules involved in reactions are small, their concentrations fluctuate significantly due to random reaction events. These concentration fluctuations cause propensity functions to also fluctuate. If the functions vary nonlinearly with concentration, the mean values of the functions have been shown to depend on the concentration fluctuations [Gomez-Uribe, et al. J. Chem. Phys. 126: 024109 (2007)]. We apply this analysis to the propensity functions to investigate various reaction network motifs and provide qualitative understanding of the fluctuation effect on bistability, stochastic focusing, and concentration detection. We have shown that this analysis allows metabolic control analysis to be applied for stochastic reaction systems by providing new versions of summation and connectivity theorems. The theorems show that sensitivities are related to one another as in metabolic control analysis. However, elasticities may not be viewed as local measures of sensitivity as they are in MCA due to noise propagation through the connected reaction pathways.
Caulobacter crescentus is an important model organism for studying the regulation of cell division cycle and cellular differentiation in prokaryotes. Caulobacter undergoes asymmetric division producing two progeny cells with identical genome but different developmental programs: the "swarmer" cell is flagellated and motile, and the "stalked" cell is sessile. Only stalked cells undergo chromosome replication and cell division. Swarmer cells must shed their flagellum and grow stalks before they can enter the replication-division cycle. Based on published experimental evidence, we propose a molecular mechanism for cell cycle control in this bacterium. Our quantitative model predicts detailed temporal dynamics of regulatory gene expression during the cell cycle and differentiation process of wild-type cells (both stalked cells and swarmer cells) as well as several mutant strains. The simulation presents a unified view of temporal regulation of protein activities during the asymmetric cell division cycle of Caulobacter. It helps to interpret phenotypes of known mutants and predict novel ones. The model can serve as a starting point for investigating the regulation of cell division and differentiation in other genera of alpha-proteobacteria, such as Brucella and Rhizobium, because recent experimental data suggest that these alpha-proteobacteria share similar genetic mechanisms for cell cycle control.
An important question in Neuroscience is how a network of neurons is affected by external stimuli. The response of the network of neurons is highly dependent on the dynamics of external stimuli as well as how they are intrinsically coupled to each other. This work is aimed at understanding how external (noisy) forces and intrinsic dynamics (coupling) effects the dynamics of a population of neurons. We assume neurons are in the oscillating regime where they are firing action potentials repetitively, and model the external forces as shared and unshared white noise. We explore the effects of weak white noise forcing (shared and unshared) on two oscillators with weak coupling. We find shared noise can induce bistability between synchronous and asynchronous states even though asynchrony is the only observable (stable) state in the absence of noise. We perform asymptotic analysis on the associated Fokker-Planck equation. To determine the particular 0th order density (\rho0) of interest, we impose the solvability condition previously outlined in (Papanicolaou 77) where we require the 1st order density to have a nontrivial solution. This requirement results in an infinite recursive relation for the coefficients of \rho0, which, a priori, is not a well-posed problem. By requiring these coefficients to be bounded, we determine conditions on the PRC and coupling function that makes this a well-posed problem. We find this \rho0 is much more computationally efficient to solve than both Monte Carlo simulations and the 2D Fokker-Planck equation; and it agrees remarkably well with Monte Carlo simulations with weak noise and weak coupling.
Boreal toads, widely distributed in the mountainous regions of the western United States, have disappeared in 80-90% of their former ranges. In this work, we model the population dynamics of boreal toads in several ways. We consider a One-Habitat, One-Sex model, a One-Habitat, Two-Sex model and a Two-Habitat, Two-Sex model that follow the developmental stages of the toad and movement between different sites. It is shown for the models that the extinction equilibrium always exists and a positive equilibrium exists under certain conditions on the population parameters. For the One-Habitat, One-Sex model, the zero equilibrium is stable if and only if the positive equilibrium does not exist. A threshold R is derived for this model such that if R < 1, the population dies out and if R > 1 the population approaches a positive equilibrium. The threshold R provides information about demographic parameters important to the preservation of the population. Apart from stabiliy, it can be shown that the One-Habitat, One-Sex model is permanent under weaker assumptions imposed on the model functions. Some numerical simulations of the three models are presented to illustrate the behavior of the population dynamics under various assumptions.
While the amount of gene expression data generated with microarrays is rapidly increasing, mathematical tools for pattern extraction from such data remain largely unchanged. In this work we introduce a new mathematical formalism that measures the extent to which genes follow a specific pattern. In particular, we focus on periodic patterns and define a family of gene periodicity measures. These measures are based on the new concept of topological persistence and are stable with respect to noise. We apply our method to the dataset of gene expressions during somite development in mouse embryo. Since this development is (approximately) periodic, genes whose expressions follow the same period are good candidates for participating in the process. We evaluate our results based on genes for which there is biological evidence of their direct involvement in somite development.
In this work we apply methods from statistical mechanics to calculate the probability that the DNA double strand opens beyond a threshold at specific locations. In particular, we evaluate partition functions and probabilities, using a nonlinear model (due to Peyrard-Bishop-Dauxois) as well as an Ising-type model (due to Frank-Kamenetskii and collaborators). Both models take into account two contributions to DNA stability, namely the hydrogen bonding between complementary bases and the stacking interaction between adjacent bases. The main purpose of these calculations is the prediction of functionally active regions along the DNA strand, based on the hypothesis that in order the bases to be read, the double strand must be open thus exposing the bases.
We compare results obtained by the two models, and also present findings on the analysis of a statistically significant number of human sequences obtained from the eukaryotic data base.
Anesthetics are believed to modulate activity of a variety of voltage-gated and ligand-gated ion channels. However, the exact biophysical mechanisms underlying anesthesia remain unclear. Recent experimental evidence suggests that the volatile anesthetic isoflurane targets the TREK and TASK two-pore potassium conductances and the persistent sodium conductance. Furthermore, recent data indicate that volatile anesthetics produce immobility, a fundamental element of anesthesia, predominantly through direct action at the level of the spinal cord.
We use mathematical modeling and experiments on the lamprey spinal cord to study the effects of isoflurane on neuronal activity. Our experimental results on the disinhibited spinal cord preparation suggest that the excitatory interneurons of the CPG are the main target of isoflurane. As the anesthetic concentration increases, ventral root activity (assumed to be representative of the activity of spinal cord excitatory interneurons) transitions from bursting to silent. However, in some animals the activity transitions directly from bursting to silent where in others it transitions from bursting to a brief period of tonic firing before falling silent.
We incorporate the anesthetic-sensitive TREK, TASK and persistent sodium conductances into a preexisting detailed biophysical model of the lamprey excitatory interneuron, and a canonical bursting model. We then perform a thorough bifurcation analysis on these models. Our results suggest that the anesthetic effects of TREK, TASK and persistent sodium currents alone are sufficient to account for the transition from bursting to silent, but are not sufficient to account for a robust transition from bursting to tonic to silent.
Many quantitative traits are affected by multiple genes. To map those quantitative trait loci (QTL) in genomic positions, we can fit data with a model that includes multiple QTL and perform a model selection analysis. For such a model selection analysis, the choice of good criteria is an open problem. A common problem in many of the methods for mapping QTL that scan positions across the genome is the determination of a threshold value for the test statistic. Many factors, such as genome size, genetic map density, informativeness of markers, and proportion of missing data, may affect the distribution of the test statistic. Permutation test has been widely used to empirically estimate the null distribution and threshold. It works well for the null of no QTL, not so well for the null of some QTL. Also it is computationally intensive. Zou et al. (2004 Genetics 168:2307-2316) proposed a model- and data-based resampling procedure using score statistics to empirically calculate the relevant threshold. It is much more efficient computationally and flexible to be adapted to complex models. We extend this procedure for multiple QTL on multiple traits to aid the model selection in the framework of multiple intervals mapping for multiple traits (MT-MIM). We have performed a simulation study to evaluate the performance of this procedure for model selection on multiple traits and the type I error has been shown to be well controlled.
Recent advances in the study of DNA indicate that three-dimensional structures other than the standard Watson-Crick double helix play roles in the mechanisms of DNA and may be responsible for genetic instability. One example of such a structure is a cruciform, which occur only when a sequence of DNA exists in inverted symmetry within the same strand. Natural stresses that occur when DNA is transcribed or replicated have been shown experimentally to drive the transitions to alternate structures by first separating the two strands of DNA. When a double stranded DNA that contains and inverted repeat (IR) becomes single stranded, the single strand can then pair with itself, forming a cruciform or hairpin structure (a.k.a. cruciform extrusion). Extrusion of cruciforms from IRs has been documented but there is a paucity of proven examples for the regulatory roles of these structures (Benham et al., 2002). It becomes useful to consider what effects their activities might have on the distributions of IRs within DNA. Many hypotheses of cruciform function, if correct, would result in IR distributions that differ from random. We use the IRFinder algorithm developed by Benson et al. to find all IRs, both those with perfect and those with approximate invert repeat symmetry, that have the attributes needed to extrude a cruciform under moderate negative superhelicity (Warburton et al., 2004). Thus, using a statistical genomics approach, we evaluate the distribution of IR sequences minimally required for cruciform formation in the Saccharomyces cerevisiae genome, finding a statistically significant enrichment of IRs that overlap or cluster at gene starts and ends, perhaps supporting a regulatory role for these structures. We also examine both inherent features of IRs (degree of imperfection, repeat length, etc) and distributional properties (whether an IR is isolated, overlap frequencies, and clustering or avoidance of regulatory regions).
Ordinary differential equations models are widely used to study IP3-mediated calcium oscillations. There, a cell is taken as a well-mixed compartment but a steep intracellular calcium concentration gradient is often observed experimentally. In addition, calcium release channels tend to be clustered at specific sites and are tightly modulated by the local calcium concentration. Using a one-dimensional spatially extended model, we investigate the effect of having spatial diffusion and discrete release sites on calcium release and oscillation. When the diffusion coefficient is large, we show that the spatial model can be reduced to an ODE model by averaging. However, calcium buffering inside the cell can reduce the diffusion coefficient value. In that case, tracking the average concentration alone is insufficient. We show that changing the value of the diffusion coefficient can significantly affect the critical IP3-level for oscillation onset as well as the frequency and amplitude of oscillation.
The intrinsic mechanisms underlying burst generation in vitro are generally well understood. However, how these mechanisms are implemented under in vivo conditions is still not clear. Pyramidal cells within the ELL have a well defined burst mechanism in vitro which is based on a somato-dendritic interaction. Surprisingly, in vivo recordings from these cells do not show any of the characteristics associated with bursting found in vitro. We propose computational model, which explain these differences.
The transient behavior of an ecosystem perturbed away from equilibrium can be described in terms of the "reactivity" and "amplification envelope" associated with the equilibrium. These properties describe the transient amplification of perturbations. They depend on the parameters of the underlying model, and sensitivity analysis predicts the effect of parameter changes on the transient response of reactive ecosystems. By doing so, it provides a quantitative framework for investigating the ecological mechanisms and processes leading to reactive dynamics.
We have developed new mathematical approaches to investigate the dependence of transient dynamics on the parameters, and will present analytical formulations for the sensitivity and elasticity of the reactivity and the amplification envelope. We will show the results of applying the methodology to a predator-prey model and a size-structured food web model. We find that perturbations can be amplified in two different ways. The predation rate may be slow but the energy transfer highly efficient, in which case the consumer can take advantage of transient increases of resources (prey-driven mechanism). Alternatively, resources may be removed quickly but poorly assimilated by the consumer, in which case the resources can grow opportunistically when the consumer biomass is below equilibrium (predator-driven mechanism). Reactivity in the food web model is highly elastic to assimilation efficiency, suggesting that the amplification of perturbations is related to the efficiency of energy transfer between trophic levels.
Many experiments have demonstrated positive effects of protist grazing on bacteria-mediated carbon decomposition. One popular hypothesis to explain this phenomenon is that nutrients recycled from the grazing of bacteria by protists facilitates carbon decomposition. Based on ecological stoichiometry, we suggest that the underlying mechanism for this phenomenon may be the cell quota difference between bacteria and protozoa. Bacteria are generally more nutrient-rich, compared to most living organisms including protozoa. This difference of cellular nutrient contents means that protozoa have to exude extra mineral nutrients to maintain cellular carbon:nutrient ratio at a constant level, i.e., strict homeostasis. Here, we develop a stoichiometry-based mathematical model to examine this hypothesis. Analytical and numerical results of our model support the hypothesis and provide us a complete view of all possible scenarios. Our results also show that under different environmental conditions, protozoa have quite different effects on carbon decomposition.
Patterns of contact in social behavior and seasonality due to environmental influences often affect the spread and persistence of diseases. Models of epidemics that assume seasonality and patterns in the contact rate include time-periodic coefficients making the systems nonautonomous. There is no general method for calculating the basic reproduction number, the threshold for disease extinction, in nonautonomous epidemic models. In some epidemic models with periodic coefficients and constant population size, the time-averaged basic reproduction number has been shown to be a threshold for disease extinction. We extend these results by showing that the time-averaged basic reproduction number is a threshold for disease extinction when the population demographics are periodic. The results are shown to hold in epidemic models with periodic demographics that include temporary immunity, quarantine, and multiple strains.
Primary growth is characterized by cell expansion facilitated by water uptake generating hydrostatic (turgor) pressure to inflate the cell, stretching the rigid cell walls. The multiple source theory of root growth is based on the proposal that phloem and protophloem transport water into the growth zone from the mature tissue higher in the root. Here protophloem water sources are used as boundary conditions in a generalized three-dimensional model of growth sustaining water potentials in primary roots. The model predicts small radial gradients in water potential, with a significant longitudinal gradient. The results improve the agreement of theory with empirical studies for water potential in the primary growth zone of roots of Zea mays. A sensitivity analysis quantifies the functional importance of small root diameter and apical phloem differentiation in permitting growth under water stress.
Syphilis was believed to disappear in U.S.A. and other developed countries after the effective treatment by penicillin years ago. However, recent resurgence of syphilis in several communities alarms a warning sign for health organizations such as CDC. Whether cycling occurs for syphilis is still on debate, however, researchers agree that dynamics of mathematical models of syphilis determines the resurgence.
We studied a few mathematical models of syphilis, incorporating the fact of partial immune protection after treatment and behavioral protection through education. In addition, like gonorrhea, treatment of patients in earlier stages may not generate any biological protection against reinfection. Our model suggests that a backward bifurcation can occur if a condition is satisfied. Our model suggests that CDC's plan of elimination of syphilis can not be achieved without enough effort on education.
Nonlinear reaction-diffusion systems which are often employed in mathematical modeling in developmental biology are usually highly stiff in both diffusion and reaction terms. Moreover, they are typically considered on multidimensional complex geometrical domains because of the complex shape of embryos. We overcome these computational challenges by combining the discontinuous Galerkin (DG) finite element methods with Strang type symmetrical operator splitting technique, on triangular meshes. This permits us to avoid directly solving a coupled nonlinear system, as necessary with the standard implicit schemes. Numerical simulations of a new biologically based system for skeletal pattern formation in the vertebrate limb, are presented to show the effects of various domain geometries on the formed patterns.
