2010 Workshop for Young Researchers in Mathematical Biology (WYRMB): Abstracts and Lecture Materials
Plenary Speakers
Future discovery and control in biology and medicine will come from the mathematical modeling of large-scale molecular biological data, such as DNA microarray data, just as Kepler discovered the laws of planetary motion by using mathematics to describe trends in astronomical data [1]. In this talk, I will demonstrate that mathematical modeling of DNA microarray data can be used to correctly predict previously unknown mechanisms that govern the activities of DNA and RNA.
First, I will describe the computational prediction of a mechanism of regulation, by developing generalizations of the matrix and tensor computations that underlie theoretical physics and using them to uncover a genome-wide pattern of correlation between DNA replication initiation and RNA expression during the cell cycle [2,3].
Second, I will describe the recent experimental verification of this computational prediction, by analyzing global expression in synchronized cultures of yeast under conditions that prevent DNA replication initiation without delaying cell cycle progression [4].
Third, I will describe the use of the singular value decomposition to uncover "asymmetric Hermite functions," a generalization of the eigenfunctions of the quantum harmonic oscillator, in genome-wide mRNA lengths distribution data [5]. These patterns might be explained by a previously undiscovered asymmetry in RNA gel electrophoresis band broadening and hint at two competing evolutionary forces that determine the lengths of gene transcripts.
Finally, I will describe ongoing work in the development of tensor algebra algorithms (as well as visual correlation tools), the integrative and comparative modeling of DNA microarray data (as well as rRNA sequence data), and the discovery of mechanisms that regulate cell division, cancer and evolution.
- Alter, PNAS 103, 16063 (2006).
- Alter & Golub, PNAS 101, 16577 (2004).
- Omberg, Golub & Alter, PNAS 104, 18371 (2007).
- Omberg, Meyerson, Kobayashi, Drury, Diffley & Alter, Nature MSB 5, 312 (2009).
- Alter & Golub, PNAS 103, 11828 (2006).
This talk consists of two parts: Pattern formation in families of microtubules under the action of kinesin and the detailed motion of kinesin along a microtubule.
Microtubules are long cylindrical structures (lengths being tens of microns and diameter approximately 25 nm) comprised of tubulin dimers, which self-assemble, 13 protofilaments being required side-to-side to form the circular cross section. In the first set of results, microtubules are represented as stiff, polar rods which are subject to diffusion in position and orientation and also subject to pair-wise interaction, mediated by kinesin molecular motors. The concentration of kinesin is represented by a parameter that feeds into the probability of an interaction occurring when two microtubules collide. The probability of an interaction also depends on the location of the collision point along the lengths of the microtubules, because kinesin accumulates at the positive end of each microtubule. With collision rules in place, Monte-Carlo simulations for large numbers of freely moving microtubules are performed, adjusting parameters for concentration of kinesin and polarity of the microtubules. From these studies, a phase diagram is produced, indicating thresholds for phase change to occur. Simulation results are compared to those from in vitro experiments.
The second part of the talk involves modeling the fine scale dynamics of a kinesin motor as it walks along a microtubule. The two heads of the kinesin molecule alternately bind and unbind to the microtubule with certain mechanisms providing a directional bias to the Brownian motion expected. One bias is the shape of the head and the shape of the binding site, along with the companion electrostatic charges. The second bias is that, utilizing ATP capture and transferal of phosphors for energy, part of the polymeric leg (neck-linker) of the bound head becomes attached towards the front of that head (the \lq\lq zipped\rq\rq\ state). The trailing head detaches from the microtubule. It then becomes subject to the biased entropic force due to the zipped state of the leading head and also preferentially (because of shape orientation) attaches in front of the currently attached head at which time ADP is released and a conformational change occurs, strengthening the binding. This motion is modeled using stochastic a differential equation. Simulations are performed with different lengths of neck-linkers and the mean speeds of progression obtained. These are compared with experimental results
The construction of flagellar motors in motile bacteria such as Salmonella is a carefully regulated genetic process. Among the structures that are built are the hook and the filament. The length of the hook is tightly controlled while the length of filaments is less so. However, if a filament is broken off it will regrow, while a broken hook will not regrow.
The question that will be addressed in this talk is how Salmonella detects and regulates the length of these structures. This is related to the more general question of how physical properties (such as size or length) can be detected by chemical signals and what those mechanisms are.
In this talk, I will present mathematical models for the regulation of hook and filament length. The model for hook length regulation is based on the hypothesis that the hook length is determined by the rate of secretion of the length regulatory molecule FliK and a cleavage reaction with the gatekeeper molecule FlhB. A stochastic model for this interaction is built and analyzed, showing excellent agreement with hook length data. The model for filament length regulation is based on the hypothesis that the growth of filaments is diffusion limited and is measured by negative feedback involving the regulatory protein FlgM. Thus, the model includes diffusion on a one-dimensional domain with a moving boundary, coupled with a negative feedback chemical network. The model shows excellent qualitative agreement with data, although there are some interesting unresolved issues related to the quantitative results.
Progress in systems biology relies on the use of mathematical and statistical models for system level studies of biological processes. Several different modeling frameworks have been used successfully, including traditional differential equations based models, a variety of stochastic models, agent-based models, and Boolean networks, to name some common ones. This talk will focus on several types of discrete models, and will describe a common mathematical approach to their comparison and analysis, which relies on computer algebra. Hence, we refer to such models as "algebraic models." The talk will present specific examples of biological systems that can be modeled and analyzed in this way.
A multispecies continuum model is developed to simulate the dynamics of cell lineages in solid tumors. The model accounts for spatiotemporally varying cell proliferation and death mediated by the heterogeneous distribution of oxygen and soluble chemical factors. Together, these regulate the rates of self-renewal and differentiation of the different cells within the lineages. As demonstrated in the talk, the feedback processes are found to play a critical role in tumor progression and the development of morphological instability.
Metabolic rate, heart rate, and lifespan depend on body size according to scaling relationships that extend over ~21 orders of magnitude and that represent diverse taxa and environments. These relationships for body mass have long been approximated by power laws, but there has been intense debate about the values of exponents (e.g., 1/4 versus 1/3). I will show for mammals that these scaling relationships exhibit systematic curvature in logarithmic space. This curvature explains why different studies find different power-law exponents. I will also show how existing optimal network theory can be modified using finite-size corrections and hydrodynamical considerations to predict curvature. I will distinguish among potential physiological mechanisms by comparing model predictions for the direction and magnitude of the curvature with results from empirical data. For the final half of the talk, I will develop modified network models to describe tumor angiogenesis and vascular structure. These new models will help to compare tumor with normal vasculature, to understand different phases (pre- and post-angiogenesis) of tumor growth, and to describe the formation of a necrotic core.
Short Talks
Evidence of rapid evolution in ecological communities has accumulated in the last thirty years, yet theory explaining the interplay between ecological and evolutionary processes on comparable time scales has not kept pace. This disparity between experiments and theory is partially due to the intractability of high dimensional systems of ordinary differential equations - even systems with two evolving species must be of at least dimension four. I will present work focusing on how the theory of slow-fast dynamical systems can be used to study predator-prey system exhibiting rapid evolution. This approach not only reduces our system back down to two dimensions, but also yields graphical techniques with predictive power about the (potentially new) qualitative dynamics a given predator-prey system can exhibit.
This work focuses on describing the scaling effects on fluid flow near the bell of the upside down jellyfish (Cassiopeia sp.) using computational fluid dynamics and live experiments. The immersed boundary method is used to simulate the bell of a jellyfish as an flexible structure coupled with a porous boundary. The porous boundary represents the oral arms which protrude over the bell, altering the flow. The effect of the oral arms on vortex formation and on volumetric flow rates are analyzed across a range of Reynolds numbers.
Mosquitoes are agile, yet surprisingly stable fliers. In the air, their long, slender legs are a visually compelling feature of their body plan. Since legs are so prominent in flight, it is unsurprising that they have long been conjectured to serve as a secondary flight control system (with the primary system being subtle modulations in wing motion). This idea has remained speculative for nearly a century.
Using high speed video to capture leg, body and wing dynamics, we develop mathematically based tracking and classification techniques to identify in-flight leg activity. By combining functional data analysis methods with a physical model-based understanding of flight force generation, we examine the importance of leg activity to in-flight maneuvering. We estimate the impact of leg activity on overall flight dynamics, examining how mosquito leg motion effects mosquito stability and turning manuevers.
The need to interpret and extract possible inferences from high-dimensional datasets has led over the past decades to the development of dimensionality reduction and data clustering techniques. Scientific and technological applications of clustering methodologies include among others biomedical imaging, data mining and bioinformatics. Current research in data clustering focuses on identifying and exploiting information on dataset geometry and on developing robust algorithms for noisy datasets. Recent approaches based on spectral graph theory have been devised to efficiently handle dataset geometries exhibiting a manifold structure, and fuzzy clustering methods have been developed that assign cluster membership probabilities to data that cannot be readily assigned to a specific cluster. In this talk, we develop a novel fuzzy spectral clustering algorithm that combines seamlessly the strengths of spectral approaches to clustering with various desirable properties of fuzzy methods. We also discuss examples of gene expression datasets for which the developed methodology outperforms other frequently used algorithms.
In the generic class of biological population models, the notion of persistence (impossibility of species extinction) is often of central interest. Persistence has an intrinsic importance in the study of animal population dynamics and in the analysis of infections spread. Moreover, in certain general settings, it has been shown that persistence precludes bistable, switch-like or oscillatory behavior of the model; this fact is of major significance to biochemical interaction networks models (e.g. metabolic pathways), where a great deal of attention has been paid recently to devising criteria that determine the possibility of such behaviors.
In this work we show that a large class of two-dimensional biological interaction network models is persistent. This result is robust, in the sense that it holds independently of the choice of parameters in the model. More precisely, we prove that any two-dimensional, endotactic, k-variable mass-action system is persistent. A k-variable mass-action system is a generalized mass-action system where rate constants are allowed to vary with time. A system is endotactic if its underlying network satisfies an easily-checked geometric property. The class of endotactic networks encompasses the well-known class of weakly reversible networks.
In the end, we use our persistence result to prove the three-dimensional case of the Global Attractor Conjecture for biochemical reaction networks, a long-standing conjecture originating in the work of Horn, Jackson and Feinberg.
This is joint work with George Craciun.
In this work, a multi-stage age-structured population dynamics model is used to determine the most efficient and least expensive procedures to reduce the population size of a vine pest. Egg and larval pesticides are often used by wine makers to contain this pest. Making the decision to apply these products or not is sometimes not easy because it depends on several variables, for example weather, price or field observations. Here, we focus on three constraints that are the action mode, the price and the efficiency of pesticides and, we determine according to them the best treatment during the pest life cycle to maximum reduce the population size of this insect. To get the optimal control on this population we study an optimization problem with constraints.
Body size has been shown to be a significant factor in shaping the structure of food webs, which are network models of the flow of energy in an ecosystem. Recent studies have shown that body size constraints can influence food web dynamics through prey preference and foraging behavior, and can thereby influence the stability of these ecosystem models. Because of its significance, we use body size as the species strategy in an evolutionary game theory approach to studying the influence of predation at individual trophic levels on evolutionarily stable strategies (ESS) in food webs.
We systematically construct small (3-5 species) food webs, and combine ecological and evolutionary dynamics using differential equation models to show how the addition of each trophic level impacts the equilibrium strategies of other species. The strategy in our model influences the intrinsic growth rate and carrying capacity of the basal (plant) species, and the interaction rates across species. We show that when a consumer is introduced, the equilibrium strategy of the basal species evolves toward a value that increases the intrinsic growth rate; however, the strength of this effect is mediated by predator species at the third trophic level. We also show how size-based prey preference can influence strategy dynamics and population sizes over long time scales. These results suggest that understanding evolution of body size is important for understanding the trophic interactions that form the basis for large-scale food web structure and function.
Michaelis-Menten kinetics are often used to describe enzyme-catalyzed reactions in biochemical models. The Michaelis-Menten approximation has been rigorously derived in the context of traditional differential equation models. In many biochemical processes, however, stochastic effects due to the random interaction of chemical species present in small numbers play an important role. To capture these effects one has to move away from differential equation models, which are continuous and deterministic, to a discrete stochastic representation. The theory underlying Gillespie stochastic simulation algorithm (SSA) provides a physical justification for stochastic models comprised of elementary reactions. Solving for the evolution of the chemical species distributions exactly is intractable for most models. Typically, the SSA is used to generate an ensemble of trajectories to establish an estimate of the distribution. The slow-scale stochastic simulation algorithm (ssSSA) uses model reduction to improve the performance of the SSA. However, the ssSSA suggests an approximation that differs from the Michaelis-Menten rate for certain enzymatic reactions. First-passage time analysis is used to examine some of these differences. The Michaelis-Menten approximation is shown to be applicable under a set of validity criteria in discrete stochastic models.
The neuronal networks of the olfactory system transduce and transform complex mixtures of odorant molecules into patterns of the neural activity representing smells. We explore two important aspects of how this process works, at the cellular and the neural circuit level, in modeling studies that produce experimental testable predictions.
- It has long been known that (in rodents) initial synaptic olfactory processing occurs in the olfactory bulb (OB) glomeruli, but the roles of various juxtaglomerular neurons is still not well understood. Recent experimental studies indicate that endogenously bursting external tufted (ET) cells -- which connect olfactory receptor neurons (ORN; OB input) to mitral cells (MC; OB output) -- play a central role in coordinating intraglomerular activity. We develop a biophysically realistic, Hodgkin-Huxley-style ET cell model that includes membrane currents known to be essential for bursting. We use specialized smooth optimization methods to study the (local) sensitivity of its functional characteristics (e.g. burst duration, interburst interval) to parameters, and statistical analyses to characterize the (global) influence of different currents.
- Odorant-evoked input to and output from the OB is temporally dynamic, and these dynamics are important in shaping odor perception. Inhalation-evoked input bursts of ORN activity occur with durations, latencies, and amplitudes that vary across glomeruli (for the same odorant) and also in individual glomeruli for different odorants, and similarly diverse activity patterns occur at the MC level. We investigate these dynamics using biophysical models of the ORN-MC and ORN-ET-MC circuits. The models’ inputs are taken from recordings of ORN calcium signals of head-fixed rats exposed to odorants and closely reproduce signals received by the real neurons. With this data-driven dynamical modeling approach, we are able to explore how the circuits’ response dynamics vary for different odorants, synaptic strengths, and intrinsic cellular parameters.
Soliton like structures called "stable droplets" are found to exist within a paradigm reaction diffusion model which can be used to describe the patterning in a number of fish species. It is straightforward to analyse this phenomenon in the case when two non-zero stable steady states are symmetric, however the asymmetric case is more challenging. We use a recently developed perturbation technique to investigate the weakly asymmetric case.
Posters
This study presents a deterministic model for theoretically assessing the transmission of bovine tuberculosis in a single on farm cattle herd. The model is analyzed to gain insights into the qualitative features of its associated equilibria. This allows for the determination of important epidemiological thresholds such as the basic reproduction number. The model is first studied when the source of infection is within the cattle herd and it is shown that it has a globally asymptotically stable disease-free equilibrium whenever the reproduction threshold is less than unity. Further, the model is studied for the case where the source of infection into the herd is also due to infection as a result of restocking of the herd with infected cattle in disease latency stages. The result shows that the disease equilibrium is globally asymptotically stable.
Recent experiments in retinal ganglion cells (RGCs) have demonstrated that pairwise maximum entropy (PME) methods can approximate observed spiking patterns to a high degree of accuracy. Through numerical simulation and analytical results on a simple class of model circuits, we show that feedforward circuits can generate a very limited departure from PME fits. The degree of departure is determined by marginal input statistics, rather than the order of input correlations. As a result, we can make quantitative predictions for the quality of PME fits in RGCs under a variety of stimulus conditions. We verify our predictions with a realistic spiking model of RGCs.
Effective regulation and appropriate activation of the immune system is traditionally understood in terms of pattern recognition mechanisms, route of infection effects and dosage dependent responses. Recently, a number of experiments have suggested that antigen kinetics may play a role in immune system decision-making as well. The mechanisms underlying this mode of immune regulation, however, are unknown. We develop a system of ordinary differential equations (ODEs) to describe the kinetics and signalling interactions that have been reported between CD4+ helper T-cells (Th). We then use this model to test whether or not the Th interaction network is capable of accurately and robustly classifying pathogens based on their relative rates of population expansion. Model analysis shows that positive and negative feedback between Th cells results in multistability of the system, and that this multistability can, indeed, generate robust classification based on pathogen growth. The Th cell interactions that we use to build our model are widely known; however independent experiments have failed to rationalise their significance. Mathematical modelling allows us to integrate together a wide array of empirical observations in order to elucidate the systems-level functioning of CD4+ helper T-cells.
This article introduces a novel model that studies the major factors jeopardizing the TB control programme in China. A previously developed two strain TB model is augmented with a class of individuals not registered under the TB control programme. The paper investigates the basic reproduction number and proves the global stability of the disease free equilibrium. The presence of three endemic equilibria is established in the model. With the help of numerical simulations a comparative study has been performed to test the validity of the model presented here to the real data available from the Ministry of Health of the People's Republic of China. Sensitivity and elasticity analysis suggest the impact of key parameters on the tuberculosis control in China.
Multicellular communities are a dominant, if not the predominant, form of bacterial growth. Growing affixed to a surface, they are termed biofilms. When growing freely suspended in aqueous environments, they are usually referred to as flocs. Flocculated growth is important in conditions as varied as bloodstream infections (where flocs can be seen under the microscope) to algal blooms (where they can be seen from low earth orbit). Understanding the distribution of floc sizes in a disperse collection of bacterial colonies is a significant experimental and theoretical challenge. One analytical approach is the application of the Smoluchowski coagulation equations, a group of PDEs that track the evolution of a particle size distribution over time.
The equations are characterized by kernels describing the result of floc collisions as well as hydrodynamic-mediated fragmentation into daughter aggregates. The post-fragmentation probability density of daughter flocs is one of the least well-understood aspects of modeling flocculation. A wide variety of functional forms have been used over the years for describing fragmentation, and few have had experimental data to aid in its construction. In this talk, we discuss the use of 3D positional data of Klebsiella pneumoniae bacterial flocs in suspension, along with the knowledge of hydrodynamic properties of a laminar flow field, to construct a probability density function of floc volumes after a fragmentation event. Computational results are provided which predict that the primary fragmentation mechanism for medium to large flocs is erosion, as opposed to the binary fragmentation mechanism (i.e. a fragmentation that results in two similarly-sized daughter flocs) that has traditionally been assumed.
The formation of the layers within the cortex of the mammalian brain has been of interest to researchers in developmental biology. Experiments on mice have found a counter-intuitive layering sequence. The layers are formed in an inside-out arrangement, with the deepest layers formed first and the outermost formed last. It is believed that the glycoprotein Reelin is important to the layering mechanism. It is observed that when this protein is absent from the mouse cortex, the layer arrangement is inverted. This resulted in a number of mechanisms being proposed for how Reelin may affect the neurons that create these layers within the cortex. However, no conclusive model has yet been determined.
Motivated by this biological problem, I will present a model for convection-dominated invasion of a spatial region by initially motile agents which are able to settle. The motion of the motile agents and their rate of settling will be affected by the local concentration of settled agents. The model can be formulated as a nonlinear partial differential equation for the time-integrated local concentration of the motile agents, from which the instantaneous density of settled agents and its long-time limit can be extracted. In the limit of zero diffusivity, the partial differential equation is of first order; for application-relevant initial and boundary-value problems, shocks arise in the time-integrated motile agent density, leading to delta-function components in the motile agent density; and there are simple solutions for a model of successive layer formation. The model is found to offer insight into biological processes involving layered growth or overlapping generations of colonization.
Although much is known about the biophysics, anatomy and physiology of basal ganglia (BG) networks, cellular and network basis of parkinsonian tremor remains a vaguely understood issue. A large body of experimental evidence supports the hypothesis, that tremor arises due to pathological interaction of potentially oscillatory cells and circuits within the loop formed by basal ganglia and thalamocortical networks. We propose a model of the BG circuitry, which helps clarify the potential mechanism of tremor genesis.
We consider a conductance-based model of subthalamo-pallidal circuits (based on the models of STN and GPe neurons) embedded into a simplified representation of thalamocortical circuit to investigate the dynamics of this loop. We modulate the connections strength of the loop to represent the modulation of basal ganglia synapses by dopamine. A number of synaptic connections in the basal ganglia are known to be presynaptically suppressed by dopamine. Thus, weaker connections in our model may correspond to healthy state or to the administration of dopaminergic medication, while stronger connections may correspond to the decreased level of dopamine in Parkinson's disease. The stronger connections (lower dopamine) favor bursting in the tremor frequency range and this result is robust with respect to what particular synaptic strength is changed. Similarly, when the feedback is completely removed, STN switches back to tonic or almost tonic firing. This happens when the feedback is interrupted either at the level of output of STN or at the level of input to pallidum and STN. However, the outcomes of these two possibilities are two different types of activity.
The proposed model supports the basal ganglia thalamocortical loop mechanism of tremor generation. Under normal conditions the behavior of this loop is not oscillatory. The dopaminergic deficit of Parkinson's disease induces tremor-related oscillations in the loop. This suggests that thalamocortical circuits (which provide a feedback to basal ganglia) are changed in Parkinsons disease in such a way as to promote the birth of oscillations. Surgical intervention or dopaminergic medication in parkinsonian patients can break or weaken pathological feedback so that basal ganglia circuits are no longer able to produce bursting. This may give the answer to the question of how tremor is suppressed.
Many experiments in the pre-Botzinger (pBC) complex of the mammalian respiratory brainstem are done in synaptic isolation. Ultimately, we would like to explore the interactions of intrinsic dynamics and coupling architecture in the context of the pBC. To make progress toward this goal, we wish to model the individual cells of the pBC. In the pBC, most individual neurons express both the persistent sodium current (INaP) and the calcium-activated nonspecific cationic current (ICAN). However, previous modeling efforts have focused on either INaP or ICAN. Here, we analyze the effects of including both INaP and ICAN within one model. Using a slow-fast decomposition, we explore the effects of varying gNaP, the conductance of INaP and gCAN, the conductance of ICAN, and explain the transition mechanisms between different regimes of tonic and bursting activities. For certain ranges of gNAP, we find that increasing gCAN from 0 to 5 yields solutions that are tonically active, then solutions that have square wave bursts, then tonic activity again, and finally solutions have bursts that exhibit depolarization block. Interestingly, we also find that the presence of INaP enhances the ability of a self-coupled cell to emit bursts featuring a period of depolarization block; previously, such bursts were only seen in the model with ICAN.
Individually, ICAN and INaP have been found to be important in the bursting rhythm in other systems, such as the Entorhinal Cortex. It is possible that both currents are present in these systems, so our analysis may be applicable to other neuronal systems as well.
My work describes how the transient dynamics of nonlinear structured population models (Integral Projection Models and Population Projection Matrices) are tied in an important way to the nonlinear fecundity function chosen. I will show that two models that have identical linear data and asymptotic population distributions can have wildly different transient behavior. This event can occur even if their (distinct) nonlinearities are both in a largely accepted class of functions that have modeled the phenomena of interest. The aforementioned outcome stresses the importance of deriving such nonlinearities with the use of biological intuition, as opposed to ad-hoc statistical analysis. This work also highlights the implicit goal of modeling population dynamics as often being too focused on one end of the time index (towards infinity), when important features of the population often happen early in time.
The initial reaction of the body to bacterial infection or severe tissue trauma is an acute inflammatory response, such response helps to annihilate threat posed by endotoxins and thus restore health. In a previous work by Roy et al, an 8-state ordinary differential equation (ODE) model of the acute inflammatory response system to endotoxin challenge was developed. Endotoxin challenges at 3, 6 and 12 mg/kg were administered to rats, and experimental data for pro-inflammatory mediators such as interleukin-6 (IL - 6) and tumor necrosis factor (TNF), as well as anti-inflammatory mediator such as interleukin-10 (IL-10) were obtained. In this work, we developed a reduced ODE model by categorizing state variables with similar behavior or functions into the same group. We employed sensitivity analysis to find parameters that were sensitive to the reduced model; the subset selection method of "SVD followed by QR factorization with column pivoting" was used for parameter identifiability analysis, we then estimated the parameters identified from subset selection.
Experimental data on endotoxin challenges at 3 and 12mg/kg were used to calibrate both models (original, and reduced models), and challenge level 6mg/kg for model prediction. Model comparison and validation were done by comparing curve fittings of the original, and the reduced models against existing experimental data, Akaike's Information Criterion (AIC) was also introduced to make quantitative comparison. Finally, results obtained from comparing both models showed that the reduced model had comparative or better performance over the original model.
In this work we develop and analyze a mathematical model describing the dynamics of infection by a virus of a host population in a freshwater environment. Our model, which consists of a system of nonlinear ordinary differential equations, includes an intrinsic quota, that is, we use a nutrient (e.g., phosphorus) as a limiting element for the host and potentially for the virus. Motivation for such a model arises from studies that raise the possibility that on one hand, viruses may be limited by phosphorus, and on the other, that they may have a role in stimulating the host to acquire the nutrient. We perform an in-depth mathematical analysis of the system including the existence and uniqueness of solutions, equilibria, asymptotic, and persistence analysis. We compare the model to experimental data, and find that biologically meaningful parameter values provide a good fit. We conclude that the mathematical model supports the hypothesized role of stored nutrient regulating the dynamics, and that coexistence of virus and host is the natural state of the system.
Throughout biology, cells and organisms use flagella and cilia to propel fluid and achieve motility. The beating of these organelles, and the corresponding ability to sense, respond to and modulate this beat is central to many processes in health and disease. While the mechanics of flagellum-fluid interaction has been the subject of extensive mathematical studies, these models have been restricted to being geometrically linear or weakly non-linear, despite the high curvatures observed physiologically. We study the effect of geometrical nonlinearity, focusing on the spermatozoon flagellum. For a wide range of physiologically relevant parameters, the non-linear model predicts that flagellar compression by the internal forces initiates an effective buckling behaviour, leading to a symmetry breaking bifurcation which causes profound and complicated changes in the waveform and swimming trajectory, as well as the breakdown of the linear theory. The emergent waveform also induces curved swimming in an otherwise symmetric system, with the swimming trajectory being sensitive to head shape - no signalling or asymmetric forces are required. We conclude that non-linear models are essential in understanding the flagellar waveform in migratory human sperm; these models will also be invaluable in understanding motile flagella and cilia in other systems.
Tumor Control Probability (TCP) is the probability that no tumor cells exist after treatment. It is used to measure the efficacy of radiation treatment and for treatment optimization. The simplest TCP models are based on statistical distributions such as Poisson distribution and Binomial distribution. More sophisticated TCP models are based on birth-death processes. These TCP formulas can be computed directly from cell population models which are differential equations for the mean expected tumor cell numbers. In this poster, I will present the some TCP models and we will use these models to study six typical treatment protocols for prostate cancer.
This paper provides analytical solutions to the generalized Fisher equation with a class of time varying diffusion coefficients. To accomplish this we use the Painlev property for partial differential equations as defined by Weiss in 1983 in The Painlev property for partial-differential equations. This was first done for Fisher's equation by Ogun and Kart in 2007; we build on this work, finding additional solutions with a weaker restriction on the trial solution. As a first step towards the stability analysis of the solutions to our generalized Fisher equation, we investigate the stability of plane wave solutions to the constant coefficient Fisher equation. This is significant to the field of Biology since Fisher's equation can be used to model the advance of an advantageous gene through a geographic region.
We present a new algorithm for the identification of bound regions from ChIP-seq experiments. Our method for identifying statistically significant peaks from read coverage is inspired by the notion of persistence in topological data analysis and provides a non-parametric approach that is robust to noise in experiments. Specifically, our method reduces the peak calling problem to the study of tree-based statistics derived from the data. We demonstrate the accuracy of our method on existing datasets, and we show that it can discover previously missed regions and can more clearly discriminate between multiple binding events. The software T-PIC (Tree shape Peak Identification for ChIP-Seq) is available at http://math.berkeley.edu/~vhower/tpic.html
Many biological systems are modeled qualitatively with discrete models. Several different modeling types have established communities in the biological sciences, including probabilistic Boolean networks, logical models, bounded petri-nets, and agent-based models. These and other discrete model types can be translated into algebraic models. Using algebraic models as a representation for discrete models allows one to apply theory from algebraic geometry and tools from computational algebra to analyze the dynamic features of such systems. Simulation has become common practice for analyzing discrete models, but most real world biological systems are far too complex to be analyzed by simulation alone. We use various abstract algebra techniques to develop algorithms and software to analyze discrete models for key dynamic features of biological relevance. All algorithms and methods are available trough a web-interface http://dvd.vbi.vt.edu/cgi-bin/git/adam.pl. Analysis of Discrete Algebraic Models (ADAM) has a 'modeler friendly' interface that allows for fast analysis of large models while requiring no understanding of the underlying mathematics or installing software. By providing a user-friendly interface to fast analysis tools, we promote the use of discrete models to model large complex systems.
The critical domain size problem determines the size of the region of habitat needed to ensure population persistence. We address the critical domain size problem for seasonally fluctuating stream environments and determine how large a reach of suitable stream habitat is needed to ensure population persistence of a stream-dwelling species. We characterize the fluctuating environments in terms of seasonal correlations between the flow, transfer rates, diffusion and settling rates, and we investigate the effect of such correlations on the critical domain size problem. We show how results for the seasonally fluctuating stream can formally connected to those for autonomous integrodifferential equations, through the appropriate weighted averaging methods.
We investigate clustering of motile and proliferative malignant glioma cells [1]. In vitro experiments in collagen gels identified a cell line that formed clusters in a region of low cell density, whereas a very similar cell line (which lacks an important mutation) did not cluster significantly. We hypothesize that the mutation affects the strength of cell-cell adhesion. We investigate this effect in a new experiment, which follows the clustering dynamics of glioma cells on a surface. We interpret our results in terms of a stochastic model and identify two mechanisms of clustering. First, there is a critical value of the strength of adhesion; above the threshold, large clusters grow from a homogeneous suspension of motile cells; below it, the system remains homogeneous, similarly to the ordinary phase separation. Second, when cells form a cluster, we have evidence that they increase their proliferation rate. We have successfully reproduced the experimental findings and found that both mechanisms are crucial for cluster formation and growth.
We examine the dynamics of spatially-extended neuronal networks with synaptic depression. First, we show how a purely excitatory network with a continuous firing rate can support oscillatory activity such as self-sustained oscillations, traveling waves, and spiral waves akin to epileptiform events observed in vitro. Then, we study standing bump solutions in a lateral inhibitory network with a Heaviside firing rate. It has been suggested that standing bumps maybe a neural substrate of working memory. Due to the piecewise smooth nature of the system, we use a novel approach to analyzing bump stability where the spectral equation form depends on the sign of perturbation. Finally, we study binocular rivalry in a competitive neural network model with synaptic depression. In particular, we consider two coupled hypercolums within primary visual cortex (V1), representing orientation selective cells responding either to left or right eye inputs, respectively. Cells with similar (dissimilar) orientation preference are more likely to excite (inhibit) one another. Intracolumnar and intercolumnar synaptic connections are modifiable by local synaptic depression. When the hypercolumns are driven by orthogonal oriented stimuli, it is possible to induce oscillations that are representative of binocular rivalry. We first analyze the occurrence of oscillations in a space-clamped version of the model using a fast-slow analysis. We then analyze the onset of oscillations in the full spatially extended system by carrying out a piecewise smooth stability analysis of single (winner-take-all) and double (fusion) bumps within the network. In regions of parameter space where double bumps are unstable and no single bumps exist, binocular rivalry can arise as a slow alternation between either population supporting a bump.
We consider a hyperbolic-parabolic system derived from the Keller-Segel model describing repulsive chemotaxis. We establish several sharp wellposedness results, large data global regularity and explicit quantitative decay rates. Several new mathematical tools are introduced to overcome the difficulties caused by weak dissipation and low frequency accumulation in low dimensions. These results (and new techniques) have applications to many important problems in mathematical biology.
Angiogenesis, the growth of new capillaries, plays an important role in cancer, ischemia, arthritis, and other diseases. During capillary growth, endothelial cells undergo a variety of behaviors such as migration, proliferation, and apoptosis. Agent-based modeling of these behaviors contributes to our understanding of angiogenesis. As with other types of models, an agent-based model allows one to experiment by changing model parameters. However, at times it is desirable to also vary the rules governing agent actions. We explore a flexible and modular multi-scale agent-based model of endothelial cell growth. Rules are promoted to the level of independent objects. This supports the rapid testing and comparison of different hypotheses of cell behavior
Integro-difference equations are used to model spatial spread of species with nonoverlapping generations. We look at a two species competition model with Ricker's growth functions in the form of integro-difference equations. We investigate spatial dynamics about how an introduced competitor spreads into a habitat pre-occupied by a resident species. We found a formula for the so called spreading speed at which the resident species retreats and the introduced species expands in space. We also obtained conditions under which the spreading speed can be characterized as the slowest speed of a class of traveling wave solutions. In addition, we conducted numerical simulations and showed that a traveling wave solution can have a complicated tail. (Joint work with Bingtuan Li)
The functions of large biological molecules are related with various kinds of large-amplitude molecular motions. With some assumptions, those motions can be analyzed by Normal Mode Analysis and Gaussian Network Model. However, despite their contributions to wide applications, the relationship between NMA and GNM requires further discussion. In this review, the author compares the Normal Mode Analysis and Gaussian Network Model and addresses their common applications in structural biology.
A multiscale model for vascular tumour growth and angiogenesis in 3D is presented. The model combines blood flow, angiogenesis, vascular remodelling, interactions between normal and tumour cells and diffusive nutrient / VEGF transport as well as cell-cycle dynamics within each cell. A key-point for drug delivery is the structure of the vascular system that differs significantly in tumours compared to normal tissue. We follow two different strategies to reproduce vascular structures in silico. The first one is poorly virtual, therein we start with a few straight initial vessels and then simulate angiogenesis that is induced by hypoxic normal or quiescent tumour cells. The resulting network ensures that the tissue segment under consideration is adequately nourished by oxygen. As a second approach we follow a hybrid strategy in which the initial vascular structure is based on a real network obtained from scanning electron microscopy from a rat brain whi ch is then completed by simulating angiogenesis.
Early clinical manifestations of Influenza A virus (IAV) infection are attributable in part to the activation of the innate immune response and inflammation, which, when very robust, may be harmful to tissue and result in organ dysfunction. This mechanism may be an important component of the morbidity associated with pandemic strains of IAV. We developed a multiscale computational model of the host-virus dynamic, extending prior work. The model introduces a dynamic inflammatory response with relevant biological compartments and expands upon existing models of innate and adaptive immune responses to IAV. The model simulates the response and dynamic regulation between virus, immune and respiratory cells, and signaling macromolecules using a compartmental system of autonomous ordinary differential equations. Simulations describe a system evolving to either recovery or death from overwhelming tissue damage, depending on initial viral load and other virus- or host-specific parameters. The model was calibrated to a dataset collected from test mice infected with A/PR/8/34. Cohorts of mice were either young adult (2-4 months) or mature adult (18-24 months) and were infected with a non-lethal or lethal viral aliquot (50 or 500 plaque forming units intranasally). Global and local sensitivity analyses are used as model reduction technique to identify highly sensitive parameters. Calibration of model parameters is obtained using a Metropolis-Hastings-based algorithm, using simulated annealing and parallel tempering to approximate the probability density function (PDF) in parameter/initial condition space.
Full marginal PDFs of parameters were compared between cohorts, yielding well-defined hypotheses as to how the host-virus response changes with age. Issues in model identification and validation are addressed. Parameter differences in key cell types demonstrate the fundamental difference in the immune responses; we form hypotheses regarding the biological system from the mathematical model. Finally, simulations were conducted comparing the effects of the inflammatory system on total viral load and survival.
The corpus luteum (CL) is an exceptionally dynamic organ which is located in the ovary. It develops from the wound space created by ovulation of a follicle, and is responsible for the progesterone production which maintains pregnancy. The growth of the bovine CL appears to be regulated primarily by the angiogenic growth factor, fibroblast growth factor-2 (FGF-2). FGF-2 stimulates capillaries to form new blood vessels which increase the supply of nutrients to the CL.
We have constructed an ordinary differential equation (ODE) model for the development of the CL in vivo. The model consists of four variables, one for FGF-2, and three for the cell types (endothelial cells, luteal cells, other cells such as pericytes) which we postulate are the main cellular components during the early CL development.
The model assumes that FGF-2 is produced by endothelial cells, and by luteal cells at a rate that depends on the vascular density. The volume occupied by each cell type is assumed to grow at a rate dependent on the vascular volume, and each cell type stops growing only when the CL volume exceeds a threshold related to the initial wound volume. This reflects the fact that the CL is to a large extent constrained to grow within that space.
For certain parameter sets, the resulting model can reproduce several characteristic features of the growth and proliferation of the cells and the production of the growth factor, in particular the surge of FGF-2 at early times, and CL growth to a steady state volume with the different cell types in appropriate proportion.
Using a combination of numerical and analytical techniques we determine the impact of varying key model parameters such as the maximal rate of the endothelial cell proliferation. If it is sufficiently small, then the vasculature remains at a low level, and, depending on the choices of other system parameters, the volume of the luteal cells or that of the other cells approaches zero. Thus, if the nutrient supply is limited, only one cell type can survive. Conversely, if the rate of endothelial cell proliferation is large enough, then the CL will be completely occupied by vessels and there will be not enough space for the luteal and the other cells to grow.
Experimental evidence indicates that an understanding of the physical organization and positioning of chromosomes in the nucleus may lead to insight into understanding gene expression. Specifically, it is believed that certain locations in the nucleus have a higher propensity of active genes. Moreover, it is believed that the relative positioning of genes within the chromosome may also play a role in gene expression. Thus, an exploration of the geometry of chromosomes should lead to a deeper understanding of how genes are activated.
To achieve this understanding, we are considering a simpler set of chromosomes than the human body: polytene chromosomes in the fruit fly, known as Drosophila. We have developed a geometrical representation of these chromosomes using the Frenet-Serret equations and techniques from imaging analysis and Fourier Transforms. In addition, a method for optimizing estimated parameters from the image is developed using methods from variational calculus.
Presented by William Erik Sherwwod
The quasi-steady state assumption (QSSA) forms the basis of a rigorous framework for modeling the dynamics of enzyme-catalyzed biochemical reactions, and it provides a solid mathematical justification for the Michaelis-Menten kinetics often used in modeling intracellular processes. A critical supposition of QSSA-based analyses is the conservation of enzyme, which is assumed neither to be added to a reaction nor to be degraded. Yet in real cells, transcription and translation of genetic elements may add to the total enzyme pool, while regulatory molecules in the cytosol may reduce it. Such considerations cast some doubt on the adequacy of the "enzymatically isolated" QSSA for analyzing biochemical processes in the context of widely varying (local) environmental conditions. We present two approaches extending the QSSA to situations involving enzyme input and degradation, and we derive conditions under which our new approximations are valid. Numerical simulations demonstrate that our open QSSA produces approximate solutions which remain quite accurate outside the provable range of validity. Our analysis also elucidates a novel connection between the behavior of non-isolated enzymatic reactions and the dynamics of a forced damped harmonic oscillator, which potentially opens a new direction for studying the dynamic behavior of more complex biochemical networks.
We study numerically and analytically networks of neural fields that are described by one dimensional reaction-diffusion equations coupled by ephaptic connections. The neural fields are coupled in several geometries like parallel coupling (as in the auditory system), feed forward coupling, branching (as in the dendritic tree), etc... Assuming symmetry/synchrony assumptions we find activity patterns in terms of traveling waves as fronts, pulses and periodic trains. In particular, we address the plausible patterns that are preferable in the spontaneous activity of the network - patterns that emerge when the input into the network is absent or small.
During mitosis chromosomes attach to dynamic microtubules in order to move to the cell equator. Kinetochores facilitate chromosome movement by maintaining a floating grip on inserted polymerizing or depolymerizing microtubules. In many vertebrates, chromosomes experience striking oscillatory movement both close to a pole and around the equator. Several proteins that can bind the microtubule lattice have been localized at kinetochores. Yet, there is no clear understanding of the force producing mechanism at the microtubule attachment sites. Here we develop a mathematical model for force generation at the microtubule-kinetochore interface. The motor is modeled using a jump-diffusion process that incorporates both biased diffusion due to microtubule lattice binding by kinetochore elements, as well as thermal ratchet forces due to microtubule polymerization against the kinetochore plate. A key model result is that kinetochore motors yield distinctively nonlinear force-velocity relations. Furthermore, for weak binding and low activation energies, we find that the motors respond with velocities that are fairly insensitive to loads and depend on the inserted polymer growth/shortening rates, in agreement with chromosome movement data. In the case of weak binding at kinetochores, the numerical results for the motor force-velocity relation as well as breaking-load calculations are in complete agreement with our derived approximate analytic solutions. Finally, we find that variations in the arrangement of kinetochore binders on the MT lattice can significantly affect motor dynamics.
A novel approach is proposed for the determination of the ensemble of structures for a protein given a set of distance bounds from NMR experiments. In this approach, similar to X-ray crystallography, a protein is assumed to have an equilibrium structure with the atoms fluctuating around their equilibrium positions. Thus, the structure determination problem can be formulated as an optimization problem, i.e., to find the equilibrium positions and maximal possible fluctuation radii for the atoms in the protein, subject to the condition that the fluctuations should be within the given distance bounds. The advantages of using this approach are the following: (1) The mathematical problem to be solved becomes more tractable, and only a single solution to an optimization problem is required, while in conventional approaches, multiple solutions are sought for a system of inequalities, which in general is intractable. (2) A single structure along with the estimates on the atomic fluctuations suffices to describe the structure and its fluctuation ranges that are underestimated in conventional approaches using a finite number of structural samples. (3) A single structural model with the fluctuation radii corresponding to the B-factors can be obtained for the full description of an NMR structure and its dynamic properties, the same as in X-ray crystallography. The formulation of the optimization problem is given. The algorithm for solving the problem is described. The test results on model proteins are presented.
The process by which honey bees choose their nesting site is a democratic one. The success of the process is dependent on whether the participating individuals form a quorum at a particular site. We present models which describe this process. If the swarm has just one site to choose from, the viability of this site is based on its quality. In this case, we show the analogy between our model and standard epidemic models, and present some analytical results. If the swarm has two sites to choose from, then there are some interesting dynamics that depend on the quality of these sites, as well as their discovery times. We present some numerical results in this case. This work is in collaboration with Andrew Nevai (University of Central Florida) and Kevin Passino (OSU).
During the last decades, progress in electroencephalographic and magneto-encephalographic studies has opened up the possibility of establishing hypotheses concerning human behavior and physiology. In this direction, novel mathematical and computational models play more and more a crucial role in allowing to gain insight into the mechanisms involved in neural processing with respect to some features exhibited in large-scale responses of the human brain.
In particular, focusing on neural correlates of auditory selective attention and habituation reflected in electroencephalographic data. We propose a mean-field modeling approach for the simulation of large–scale neural correlates of selective attention and habituation neural correlates. We show that our simulations have confirmed experimental data, while also been in agreement with leading bottom-up top down theories of human brain processing.
Dividing cells maintain a stable size from one generation to the next. This suggests that they contain homeostatic mechanisms in which the division cycle is triggered when a particular size is attained. However, the biochemical details of these mechanisms remain unknown. Sensing mechanisms appear restricted to monitoring concentration changes, so how can such changes indicate cell volume? Volume and concentration are different types of quantities; the former is sensitive to changes in scale while the latter is not. This issue has been discussed and possibilities have been proposed, most of which involve measuring the time required for a cellular component to reach a critical concentration above which mitosis is triggered. Previous mathematical models have described the control of the G2/M transition. Although some describe the main cell cycle proteins in detail, none includes a specific mechanism for measuring cell size. Here we propose a minimal complexity 1D reaction-diffusion-convection mathematical model for a cell-size checkpoint based on the recently proposed mechanism. By minimal, we mean the fewest assumptions, reactions and components required to exhibit checkpoint behavior as arising from the spatial cellular dynamics of Pom1 and Cdr2 proteins during interphase. The model combines known chemical features of Pom1 and Cdr2 with the known dynamics of microtubules implemented within a growing domain framework. The model reproduces phenotypes of a mutant fission yeast strain as well as the effects of two drugs. Our simulations demonstrate that the proposed checkpoint mechanism is feasible from a quantitative perspective.
A small number of BMP molecules are involved in the dorsal pattern formation process of Drosophila, which could lead to large stochastic fluctuations. However, type IV collagen protein can locate BMP molecules and Sog around the receptor surface, which facilitates the formation of Sog/BMP and bound receptors, attenuating noise. Similar to collagen, SBP, which is produced by a positive feedback loop, can also enhance the formation of bound receptors but can remove BMP molecules through internalization. We will show the role of collagen and SBP in the temporal evolution of noise in the pattern formation process numerically and analytically.
