2007 Workshop for Young Researchers in Mathematical Biology (WYRMB): Abstracts and Lecture Materials
Plenary Speakers
Complex fluid-structure interactions are central to mammalian fertilization. Motile spermatozoa, muscular contractions of the uterus and oviduct, as well as ciliary beating generate forces that drive fluid motion. At the same time, the dynamic shapes of these biostructures are determined by the fluid mechanics. In this talk we will give an overview of the classical work in fluid dynamics that has been applied to reproduction. We will also present recent computational models that promise to provide insight into these complex, coupled dynamical systems.
Even though much progress has been made in main stream experimental cancer research at the molecular level, traditional methodologies alone are insufficient to resolve many important conceptual issues in cancer biology. For example, for the most part, it is still unknown how cancer originates, what drives its progression, and how treatment failure can be prevented. In this talk, I will describe novel mathematical tools which help obtain new insights into these processes. I will also show how the mathematical insights are combined with experimental studies through collaborations with cancer biologists. The main idea is to study cancer as an evolutionary dynamical system on a selection-mutation network. I will discuss the following topics: Stem cells and tissue architecture; Cancer and aging, and Drug resistance in cancer.
The question of correlations between relatives dates back to Francis Galton and Karl Pearson. Focusing on continuous traits, Galton was a pioneer of regression and correlation analysis. He failed to rediscover Mendel's discrete theory of inheritance. Galton's approach was championed by Pearson, who founded biometrical genetics. In 1918 R.A. Fisher was able to reconcile the biometrical and Mendelian approaches to genetics.
In this talk I will survey recent extensions to Fisher's calculation of genetic variances and covariances. The first extension is to X-linked traits in the presence of X inactivation. The second extension is to relatives derived from inbred strains. All members of an inbred strain of animals are genetically identical and completely homozygous. When animals from multiple inbred strains are mated, their children and later descendants display complicated covariances and variances. To calculate these quantities requires new combinatorial constructs. Both of these extensions serve to control background polygenic variation in QTL (quantitative trait loci) mapping. The third extension is to the estimation of the standard combinational coefficients for outbred populations from DNA chip data. These estimates offer substitutes for reliable pedigree structures in human QTL mapping.
The brain is a complex organ composed largely of neurons, glial cells, and blood vessels. In spite of an enormous literature, experimental and theoretical, on the brain, we do not have a good understanding of how it functions on a gross mechanistic level. In general, the brain maintains a homeostatic state with relatively small ion concentration changes, the major ions being sodium, potassium, and chloride, and a very important ion, calcium. I believe we can learn a lot about the brain by studying extreme phenomena, and one such phenomenon in the brain is called cortical spreading depression (SD for short). SD was discovered over 60 years ago by A.A.P. Leao, a Brazilian physiologist doing his PhD research on epilepsy at Harvard University. It is characterized by nonlinear chemical waves that propagate at speeds of 1-15 mm/min in the cortex of different brain structures in various animals, including humans, and results in massive changes in ion concentrations. For example, normal extracellular potassium is normally about 2-5 mM, but during an SD wave, it can reach as high as 30-40 mM. In humans, SD is associated with migraine with aura. To date we do not have a good explanation of how SD occurs, although a number of mechanisms have been hypothesized to be important for SD wave propagation. In this talk, I will review some of the characteristics of SD wave propagation, and examine some of the mechanisms that are believed to be important during SD, including ion diffusion, membrane ionic currents, osmotic effects, the spatial buffer mechanism, neurotransmitter substances, gap junctions, metabolic pumps, and synaptic connections. In this talk, continuum models of SD, consisting of coupled nonlinear diffusion equations for the ion concentrations, will be described.
In this lecture I will present the fundamental principles of a new field called viral dynamics. I will illustrate those principles using examples from hepatitis C virus infection. I will show how, using simple dynamic models combined with experimental data, one is able to uncover basic features of the biology of viral infection and access the efficacy of different treatment protocols. Further, I will show how the methods allow one to answer basic questions about how different antiviral agents act, e.g., do they block viral infection or do they interfere with viral replication after a cell is infected? Lastly, I will discuss how using kinetic information has changed the way HCV infected patients are treated.
Complex networks of interacting proteins control the physiological properties of a cell (metabolism, reproduction, motility, signaling, etc.). Intuitive reasoning about these networks is often sufficient to guide the next experiment, and a cartoon drawing of a network can be useful in codifying the results of hundreds of observations. But what tools are available for understanding the rich dynamical repertoire of such control systems? Why does a control system behave the way it does? What other behaviors are possible? How do these behaviors depend on the genetic and biochemical parameters of the system (gene dosage, enzymatic rate constants, equilibrium binding constants, etc)? Using basic principles of biochemical kinetics, we convert network diagrams into sets of ordinary differential equations and then explore their solutions by analytical and computational methods. We illustrate this approach with a mathematical model of cell cycle transitions in eukaryotes, based on a molecular network controlling the activity of cyclin-dependent kinase (Cdk). In this model, arrest points in the cell cycle correspond to stable steady states of the control system. As biochemical parameters of the control system change, these arrest points are imposed or lifted by transitions called "bifurcations." During normal growth and division, cell size is the critical parameter that drives progression from G1 to S/G2 to M phase and back to G1. Simple diagrams, which correlate Cdk activity with cell growth, give a new way of thinking about cell cycle control, particularly the role of checkpoint pathways in arresting the cycle. The method is generally applicable to any complex gene-protein network that regulates some behavior of a living cell.
References:
- Tyson, Chen & Novak (2001) Nature Rev. Mol. Cell Biol. 2:908-916.
- Tyson, Csikasz-Nagy & Novak (2002) BioEssays 24:1095-1109.
- Tyson, Chen & Novak (2003) Curr. Opin. Cell Biol. 15:221-231.
- Csikasz-Nagy et al. (2006) Biophys. J. 90:4361-4379.
Poster Presentations
The exponential growth of sequence databases poses a ma jor challenge to bioinformatics tools for querying alignment and annotation databases. There is a pressing need for methods for finding overlapping sequence intervals that are highly scalable to database size, query interval size, result size, and construction / updating of the interval database. We have developed a new interval database representation, the Nested Containment List (NCList), whose query time is O(n + log N ), where N is the database size and n is the size of the result set. In all cases tested this query algorithm is 5 - 500 fold faster than other indexing methods, such as MySQL multi-column indexing, MySQL binning, and R-Tree indexing. We provide performance comparisons both in simulated datasets and real-world genome alignment databases, across a wide range of database sizes and query interval widths. We also present an in-place NCList construction algorithm that yields database construction times that are approximately 100-fold faster than other methods available. The NCList data structure appears to provide a useful foundation for highly scalable interval database applications. NCList data structure is part of Pygr, a bioinformatics graph database library, available at http://sourceforge.net/projects/pygr.
The screwworm fly is a parasite that causes myiasis (larval infestations in tissues) in wounded mammals. It was eradicated from the United States, Mexico and parts of Central America by the sterile insect release method (SIRM). A permanent sterile barrier zone is now maintained in Panama to prevent renewed invasion into eradicated territory. We have modeled SIRM control of the screwworm fly with a system of reaction-diffusion equations. Our results suggest that the barrier zone could be shortened substantially, reducing costs without risk of screwworm reinvasion.
Identification of detailed features of neuronal systems is an important challenge in the biosciences today. Cilia are long thin structures that extend from the olfactory receptor neurons into the nasal mucus. Transduction of an odor into an electrical signal occurs in the membranes of the cilia. The Cl(Ca) channels which reside in the ciliary membrane are activated by calcium and allow a depolarizing efflux of chloride and are thought to amplify the electrical signal to the brain.
A mathematical model consisting of primarily partial differential equations is developed to model experiments, one involving the interplay between the CNG amd Cl(Ca) channels and other involving the diffusion of calcium into cilia and the resulting electrical activity. The unknowns in the problem are the concentration of calcium, a buffer, the membrane potential and, the quantity of most interest in this work, the distribution of Cl(Ca) ion channels along the length of a cilium. A simple numerical method is derived that can be used to obtain estimates of the spatial distribution of Cl(Ca) ion channels along the length of a cilium.
Spreading Depression (SD) was first observed by Leao in 1944. It is a slowly traveling wave phenomenon elicited in the cortex of various brain structures in many different animals. It is characterized by depression of electroencephalographic activity. The wave is accompanied by increased blood flow and is followed by a period of vasodilation. The wave propagates with a speed of 1-15 mm/min. SD is believed to be the underlying physiological process that causes migraine with aura. The various factors affecting SD are ionic currents, pumps, neurotransmitters, spatial buffering, and osmosis. We consider the previous model of SD by Shapiro (2000) and try to simplify his complex model. We then will introduce additional factors that are believed to affect SD and build a new model. The model consists of nonlinear partial differential equations which are to be solved numerically.
The spread of infectious disease through communities fundamentally depends on the underlying patterns of contacts between individuals. Generally, the more contacts one has, the more vulnerable one is to infection during an epidemic. Thus, outbreaks disproportionately impact the most highly connected demographics. Epidemics can then lead to sparser networks, through immunization or removal of individuals, which are more resistant to future transmission of a given disease. Using several classes of contact networks Poisson, scale-free, and small-world we characterize the structural evolution of a network due to an epidemic in terms of frailty the degree to which highly connected individuals are more vulnerable to infection and interference the extent to which the epidemic cuts off connectivity among the susceptible population that remains following an epidemic. The evolution of the susceptible network over the course of an epidemic differs among cl! asses of networks; frailty, relative to interference, accounts for an increasing component of network evolution on networks with greater variance in contacts. The result is that immunization due to prior epidemics can provide greater community protection than random vaccination on networks with heterogeneous contact patterns, while the reverse is true for highly structured populations.
Carcinogenesis is the process by which healthy cells start growing out of control. Through a series of genetic mutations these cells become cancer cells and become capable metastasis and of invading other tissues. Although mutations are random at the genetic level, at the phenotypic level most cancers evolve in a 'characteristic' if not deterministic way. We will show a Cellular Automata based model that can be used to explain this different behaviour as well as the influence of the microenvironment in the process of cancer evolution.
A two drug interaction is usually predicted by individual drug pharmacokinetic. An innovative drug-drug interaction method based on a three-level hierarchical Bayesian meta-analysis model is developed that includes a Monte Carlo Markov chain pharmacokinetic parameter estimation procedure. Underlying the parameter estimation procedure is a fast integration method of the stiff pharmacokinetic equations. In this poster, we report the results of the fast integration method and comparisons with standard stiff solvers.
Klebsiella pneumoniae and Staphylococcus epidermidis are the most common causes of intravascular catheter infections. Given time, infected devices in the bloodstream become a source of a bloodborne plume of mediators, bacteria, and bacterial and host matrix. The dislodged material can actually leave the catheter surface at nearly half a meter per second, either coming to rest in a microvascular debris field in the lung or passing through into the arterial circulation. Our current model for the dynamics of the size-structured population of aggregates in a flowing system is based on the Smoluchowski coagulation equations.
In this poster, we explain the source of the apparent "waves" in time-series distributions and their relationship to multiple bacterial phenotypes. We also rigorously derive an alternative form for a fragmentation kernel for bacterial aggregates in laminar flow. We use these results to predict the persistence of capillary-sized aggregates in circulation.
T cells are an indispensible part of the adaptive immune system. The body maintains a pool of naïve T cells a defense against infection throughout life. T cell development is an ongoing process occuring in the thymus, an organ located near the heart. We explore immature T cell differentiation and migration dynamics using systems of ordinary and partial differential equations. Our current focus is on the seeding of the thymus and the outward migration of thymocytes into the outer region of the thymus. Some insights into the constant output of mature T cells by periodic seeding of the thymus are gained.
A major challenge of Systems Biology is to reconstruct biological networks directly from data, particularly from large-scale 'omics' studies. Several methods have been proposed to reverse engineer biological networks from data using a variety of theoretical frameworks, such as statistical analyses, machine learning, chemical kinetics, metabolic control analysis, or algebra, among others. These algorithms are also very diverse in terms of the type, and amount of data that they need as input, many of them requiring completely different experiments from the others. In general, the proponents of the methods demonstrate their performance on available experimental data or on simulated data, and comparisons with other methods are usually limited. Nevertheless it is hard to objectively compare the perfomance of each method against the others due to several issues; a) many methods require different types of experiments to be performed, and in general each of them has focused on different networks; b) simulated data is often produced in ways that favor the algorithm proposed; c) simulated data is often produced with an in silico network that is very simplistic and unchallenging.
It is the aim of this work to perform this comparison in a systematic way, by determining how well each of the methods perform on a given data set. To achieve this, a synhetic biological system will be used. The use of mathematical models of biological systems allows one to simulate all sorts of conditions. The model used in this work was developed by one of us and is intended to portray all levels of biological organization (genes, proteins and metabolites), being comprised of 59 variables. The model will allow us to obtain data of any type (time-course or steady state) and will enable us to perform any kind of experiment that is required (environmental perturbations, variable perturbations and knock-out experiments), thus making it possible to assess how close or how far each method gets to a good description of the network.
Quantifying the impacts of private owned land on local and regional ecological processes is important for protecting essential ecological function of the biosphere. In this study, we studied the monthly Normalized Difference Vegetation Index (NDVI) dynamics by SPOT VEGETATION from April, 1998 to December, 2004 on three adjacent forested areas with same forest type but different proportions of private land at the Bankhead National Forest of Alabama, USA. Forested area a, b, and c is covered by 6%, 20% and 55% of private land, respectively. Higher proportion of private land (e.g., 55%) among public land can decrease annual mean NDVI values, coefficient of variance, seasonal maximum NDVI values and absolute value of rate of NDVI increase/ rate of decrease, but increase seasonal minimum NDVI. There exists spatial synchrony of NDVI dynamics between the area a and b, but it disappears between the area a and c. Information entropy at multiple temporal scales, which is related to the temporal complexity of NDVI dynamics, shows that for each single area the dynamics of information entropy across temporal scales is quite similar and there are two domains and one transition zone for scaling relationship, but it changes dramatically below the temporal scale of 4 months if area with more private land is incorporated, such as the area b and c. Our results indicate that a higher proportion of private forest land could affect the regional forest NDVI dynamics in complex and ecologically significant ways. Maintaining proper proportion of private forest land at regional level could optimize ecological functions of forests.
Yeast cells respond to spatial gradients of mating pheromones by polarizing and projecting up the gradient toward the source. It is thought that they employ a spatial sensing mechanism in which the cell compares the concentration of pheromone at different points on the cell surface and determines the maximum point, where the projection forms. Here we constructed the first spatial mathematical model of the yeast pheromone response that describes the dynamics of the heterotrimeric and Cdc42p G-protein cycles, which are linked in a cascade. Two key performance objectives of this system are (1) amplification - converting a shallow external gradient of ligand to a steep internal gradient of protein components and (2) tracking - following changes in gradient direction. We used simulations to investigate amplification mechanisms that allow tracking. We identified specific strategies for regulating the spatial dynamics of the protein components (i.e. their changing location in the cell) that would enable the cell to achieve both objectives.
Different species of malaria parasites have long been recognised to have preferences for different age-classes of cells and it has been postulated that this might account for some of the differences in infection kinetics between the parasite species. We formulated an age structured model of red blood cell infection kinetics which incorporated a preference of the parasite to invade cells of a pre-defined age. We then applied this model to data collected from P.berghei infected BALB/C mice and were able to explain two commonly observed, yet unexplained, phenomena. It has been noted that malaria infection results in a large reduction in the number of circulating reticulocytes compared to the number observed in an equivalent chemically induced anaemia. However in our data a higher percentage of infected reticulocytes than erythrocytes was also observed. By applying our mathematical model we were able to quantify the preference of P.berghei for reticulocytes at approximately 150 fold over erythrocytes, and show that the reduction in circulating reticulocytes during P.berghei infection was entirely accounted for by preferential parasitisation and lysis. Based on our analysis we hypothesize that the bone marrow production of reticulocytes is not suppressed during a malarial episode, rather, their numbers in circulation are simply decreased by their preferential parasitisation and hence destruction. Our model is also able to explain the dip in parasitaemia after approximately 10 days of infection which is seen in almost all P.berghei infected mice, again by incorporating the strong preference of the P.berghei parasite for reticuloytes. In light of the unexpectedly strong preference uncovered of P.berghei for reticulocytes we examined more generally the effect of preferential infection of reticulocytes by a malaria parasite. We were able to show that this preference may be far higher in all species than previously thought.
The vertebrate body plan is characterized by a segmented organization, which is particularly obvious at the level of the vertebral column by the repetition of the vertebra. This basic organization is first established during embryogenesis through the periodic formation of embryonic segments called somites. This process is associated with a molecular oscillator - the segmentation clock - which has been characterized in different vertebrate species ranging from fish to mouse: this clock drives the periodic expression of so called cyclic genes in the tissue precursor of the somites, the presomitic mesoderm (PSM). So far, a handful of cyclic genes have been identified and are mostly linked to the Notch pathway.
We have undertaken a microarray approach to search for new cyclic genes, by analyzing time series of mouse PSM during one clock oscillation cycle. Due to the limiting amount of starting RNA (50ng), our protocol includes two rounds of amplification, controlled for its reproducibility and quality. We performed a mathematical (Lomb-Scargle) analysis to identify periodic expression profiles. As a validation of the method, known cyclic genes ranked among the most significant candidate cyclic genes. By hierarchically clustering the expression profiles, we identified a first cluster containing new cyclic genes related to Notch but also unexpectedly to FGF signaling. The second one contained almost exclusively new cyclic genes related to Wnt signaling, cycling in antiphase. The most significant candidate cyclic genes were validated by in situ hybridization.
As a conclusion, this first quantitative and systematic analysis of the segmentation clock system by microarrays has allowed the identification of a complex oscillating network of signaling genes underlying the mouse segmentation clock.
Localized calcium elevations known as calcium puffs and sparks are cellular signals that arise from the cooperative activity of clusters of inositol 1,4,5-trisphosphate receptors and ryanodine receptors clustered at calcium release sites on the surface of the endoplasmic reticulum or sarcoplasmic reticulum. When Markov chain models of these intracellular calcium-regulated calcium channels are coupled via a mathematical representation of local calcium signaling, simulated calcium release sites may exhibit the phenomenon of "stochastic calcium excitability" where the inositol 1,4,5-trisphosphate receptors (IP3Rs) or ryanodine receptors (RyRs) open and close in a concerted fashion. Interestingly, some single channel models of calcium-regulated intracellular channels that include calcium inactivition are not particularly sensitive to the density of channels at simulated calcium release sites, so long as the requirement for inter-channel communication is satisfied, while other single channel models that do not include calcium inactivation will open and close in a synchronous fashion only when the channel density is in a prescribed range. This observation led us to hypothesize that single channel models that include calcium inactivation would be less sensitive to the details of release site ultrastructure than single channel models that lack a slow calcium inactivation process. To determine if this was the case, we simulated model calcium release sites composed of instantaneously-coupled calcium-regulated calcium-channels whose spatial locations were chosen from a uniform distribution on a disc of specified radius and corresponding simulations with channels arranged on hexagonal lattices with the same average inter-channel distance. Analysis of puff/spark statistics such as the index of dispersion of the steady-state fraction of open channels in the release site confirmed our hypothesis that the synchronous gating of clusters of calcium-regulated calcium channels is less sensitive to the spatial organization of release sites when the single channel model includes a slow calcium inactivation process. We also investigate the validity of several different mean-field approximations to calcium release site dynamics that do not explicitly account for the details of release site ultrastructure. The most successful approach maintains a distinction between the calcium nanodomains - each channel's substantial influence on its own stochastic gating - and the collective contribution of elevated calcium concentration from neighboring channels.
Bats are natural reservoirs of rabies. We address the maintenance of the disease in bat colonies by developing individual and population models that generate indicators of risk of rabies to bats, that provide dynamic estimates of effects of rabies on population densities, and that suggest consequences of viral exposures and infections in bats relative to physiological and ecological characteristics of bats in different habitats. We present individual models (within host) for the immune responses to a rabies virus challenge, an immunotypic disease model that describes the evolution of the disease and a disease demographics model, which is structured by immunotypic response governed by immune system efficiency. Model simulations are consistent with available data, characterized by relatively low prevalence of the virus in colonies and much higher prevalence of rabies virus-neutralizing antibodies. Under model conditions, there is a robust non-clinical state that can be attained by the exposed individual that allows persistence of the disease in the population.
The enormous accumulation of experimental data representing the activities of the living cell has triggered an increasing interest in the reverse engineering of biological networks from data. In particular, construction of discrete models for reverse engineering of biological networks is receiving attention, with the goal of providing a coarse-grained description of such networks. The poster considers the modeling framework of polynomial dynamical systems over finite fields constructed from experimental data. We present and propose a solution to a problem inherent in this modeling method: the selection of a particular polynomial model from the set of all models that fit the data. This issue is closely related to the problem of finding all polynomials that vanish on the data, the so-called ideal of points. Polynomial ideals can be represented through a special finite generating set, known as Groebner basis, that possesses some desirable properties. For a given ideal, however, the Groebner basis may not be unique since its computation depends on the choice of leading terms for the multivariate polynomials in the ideal. The correspondence between data points and uniqueness of Groebner bases is studied. More specifically, an algorithm is developed for finding all minimal sets of points that, added to the given set, have a corresponding ideal of points with a unique Groebner basis.
We consider an optimal fishery harvesting problem using a semilinear elliptic PDE model, which has logistic growth and the harvest depends on the location of the fish. We consider two objective functionals: maximizing the yield and minimizing the cost or the variation in the fishing effort (control). Existence, necessary conditions and uniqueness for the optimal harvesting control for both cases are established. The optimal control when minimizing the variation is characterized by a variational inequality instead of the usual algebraic characterization, which involves the solutions of an optimality system of nonlinear elliptic partial differential equations. Numerical examples are given to illustrate the results for both cases.
Globus pallidal (GPe) nucleus neurons play an important role in the motor control mechanism of basal ganglia circuitry under Parkinsonian conditions. We study phaselocking properties of the endogenously oscillating (period T0) GPe neurons to external excitatory synaptic input (with period Tf). Analytical conditions for 1:1 synchronization of the neuron with the external period are derived both for T0 >Tf and T0< Tf based on approximating the neuron's phase response curves with piece-wise linear temporal profiles. Simulation results with the actual model are compared. Other phaselocking regimes for periodic synaptic alpha function input are studied and explained in terms of the gating variables. For stronger synaptic input, the neuron can display a pause or cessation of spiking for time spans that are several times T0. The effect of such pauses on the output firing rate variability are studied.
Work done in collaboration with Charles J. Wilson.
When local treatment of breast cancer moved away from mastectomy to the less radical lumpectomy, it was realised that radiotherapy was essential to reduce local recurrence rates. We have developed a mathematical model of breast tumour initiation and development combined with surgery and radiotherapy simulations. Our simulations show that indeed all stray tumour cells in the tumour bed are likely to be killed by post-surgery irradiation. However, this radiation may not only eradicate tumour cells, but also harm non-tumour cells in the tissue. Using our model we present results that show that irradiation induced mutations in normal cells may give rise to a new tumour. By incoporating this result we will also examine a range of irradition stratageis that give the best chance for treatment success. Additionally, we discuss the hazards of low irradiation screening protocols that are aimed at detecting breast cancer early.
Cell colonies of bacteria, tumour cells and fungi exhibit, under nutrient limited growth conditions, complex branched growth patterns. In order to investigate this phenomenon we present a simple hybrid cellular automaton model of cell colony growth. In the model the growth of the colony is limited by a nutrient that is consumed by the cells and which inhibits cell division if it falls below a certain threshold. Using this model we have investigated how the nutrient consumption rate of the cells affects the growth dynamics of the colony. We found that for low consumption rates the colony takes on a Eden-like morphology, while for higher consumption rates the morphology of the colony is branched with a fractal geometry. These findings are in agreement with previous results, but the simplicity of the model presented here allows for a linear stability analysis of the system. By observing that the local growth of the colony is proportional to the flux of the nutrient we derive a dispersion relation for the growth of the colony interface. This dispersion relation shows that the stability of the growth depends on how far the nutrient penetrates into the colony. For low nutrient consumption rates the penetration distance is large, which stabilises the growth, while for high consumption rates the penetration distance is small, which leads to unstable branched growth. When the penetration distance vanishes the dispersion relation is reduced to the one describing Laplacian growth without ultra-violet regularisation. The dispersion relation was verified by measuring how the average branch width depends on the consumption rate of the cells and shows good agreement between theory and simulations.
We propose a novel biologically constrained three-phase model of the brain microstructure. Designing a realistic model is tantamount to a packing problem and, for this reason, a number of techniques from the theory of random heterogeneous materials can be brought to bear on this problem. Our analysis strongly suggests that previously developed two-phase models in which brain cells are packed in the extracellular space are insufficient representations of the brain microstructure. These models either do not preserve realistic geometric and topological features of brain tissue, or preserve these properties while overestimating the brain's effective diffusivity, an average measure of the underlying microstructure. For this reason, we propose that the brain is best modeled as a three-phase random heterogeneous material, where the extracellular matrix is the key third phase that is unaccounted for in conventional two-phase models. Using accurate first-passage-time techniques, we demonstrate that at the appropriate extracellular matrix concentration, the proposed biologically constrained three-phase model conforms to the diffusion properties of the brain microstructure.
During electrical excitation-contraction coupling in cardiac myocytes a small amount of calcium (Ca2+) enters the myoplasm via L-type Ca2+ channels (dihydropyridine receptors, DHPR) into a restricted space between the sarcolemma and sarcoplasmic reticulum (SR) membrane. In this `diadic subspace' the small amount of `trigger' Ca2+ causes the release of a much larger amount of Ca2+ from the SR via Ca2+-regulated ryanodine receptors (RyRs) resulting in a local increase in Ca2+ referred to as a Ca2+ spark. It is known that RyRs form near crystalline lattices on the SR membrane [Yin et al., 2005] and neighboring channels make physical contact with each other. This structural evidence in addition to planar lipid bilayer experiments showing two or more RyRs may open and close in a concerted manner, a phenomenon requiring the presence of the regulatory protein FKBP 12.6 [Marx et al., 2001], suggests channel-to-channel allosteric interactions may contribute significantly to the dynamics of spark formation and collapse.
Here the effects allosteric interactions have on spark dynamics are investigated by simulating a Markov chain model representing a Ca2+ release site composed of a cluster of minimal two-state Ca2+-activated RyRs organized on a square lattice. Model RyRs experience instantaneous Ca coupling as well as coupling via nearest-neighbor allosteric interactions.
Allosteric interactions can be incorporated as to stabilize coupled channels in the closed state only, the open state only, or both the closed and open state in a balanced fashion. Introducing allosteric coupling in any of these three ways enhances spark-like excitability in Monte Carlo simulations. However, increasing the strength of allosteric interactions to stabilize coupled channels in the closed state reduces the sensitivity of sparks to the strength of Ca2+ coupling, while increasing the strength of allosteric coupling between open channels increases spark sensitivity to the strength of Ca2+ coupling. The effects of randomly removing an increasing fraction of channel-to-channel allosteric couplings on spark duration, inter-spark interval, and spark frequency is investigated. Notably, if allosteric interactions stabilize the closed state only, this procedure increases spark duration, decreases inter-spark interval, and at first increases but ultimately decreases spark frequency.
A mean-field approximation applicable to the dynamics of a RyR cluster coupled via both Ca2+ release and nearest-neighbor allosteric interactions is derived and validated against Monte Carlo simulations. Importantly, the relationship between lumped rates in this reduced model and microscopic parameters such as the strength of cooperative allosteric coupling are preserved.
Many activities of cells are controlled by cell surface receptors, which respond to ligands by triggering intracellular signaling reactions. Improved understanding of receptor signaling has several applications, such as the rational design of drugs. The process of receptor signaling involves highly connected networks of interacting components. Understanding the behavior of these networks requires the development of mathematical and computational modeling.
Here, an agent-based simulator is implemented to study diverse molecular interactions in complex signaling pathways with spatial resolution and single molecule detail. It permits stochastic modeling of protein clustering, protein motion and biochemical reactions within an idealized cellular geometry. Components diffuse and react with nearby particles in accord with chemical rate reactions. Reactions occurring between two molecules depend on the specific types of the molecules and their physical states. The modularized design confers flexibility.
The model is presently being used to determine the impact of spatial proximity and microdomain organization upon EGFR/ErbB receptor signal transduction efficiency and outcome. Simulation results are compared with experimental data to provide insight into the details of cell membrane organization.
A numerical scheme for diffusion-advection of solutes in the solute-fluid-structure interaction is proposed. When concentrations inside and outside of the boundary are different, the moving boundary drives the solutes to keep the law of conservation. For numerical accuracy and efficiency, the fine meshing along the boundary is used. Implicit Stokes equation is solved adaptively with cell-centered discretization and approximate projection method. Diffusion-advection equations are solved with FAC and algebraic multigrid methods.
As an application of this mathematical framework, the modeing and simulation of the dendritic spine motility based on the acto-myosin network are done. The PDE version of the calcium induced calcium release, the contraction dependent on the Ca++ and IP3 concentration and the spine geometry are considered. In the case of several ions staying together, the electrical potential and charge density distribution are checked to satisfy electrical neutrality and thermodynamic relation.
Computational techniques within the population density function (PDF) framework have provided time-saving alternatives to classical Monte Carlo simulations of neural network activity. Efficiency of the PDF method is lost as the underlying neuron model is made more realistic and the number of state variables increases. In a detailed theoretical and computational study, we elucidate strengths and weaknesses of dimension reduction by a particular moment closure method (Cai et al., 2004, 2006) as applied to integrate-and-fire neurons that receive excitatory synaptic input only. When the unitary postsynaptic conductance event has a single-exponential time course, the evolution equation for the PDF is a partial differential integral equation in 2 state variables, voltage and excitatory conductance. In the moment closure method, one approximates the conditional centered moment of excitatory conductance given voltage by the corresponding unconditioned moment. The result is a system of coupled (1-D) partial differential equations in voltage. Moment closure at k=2 works well, and at k=3 works even better, in the regime of high dynamically varying synaptic input rates. Both closures break down at lower synaptic input rates. Phase plane analysis of the k=2 problem with typical parameters proves, and reveals why, no steady-state solutions exist below a synaptic input rate that gives a firing rate of 59 1/s in the full 2-D problem. Closure at k=3 fails for similar reasons. Low firing rate solutions can only be obtained with parameters for the amplitude and/or kinetics of the unitary postsynaptic conductance event that are on the edge of the physiological range. We conclude that this dimension reduction method gives ill-posed problems for a wide range of physiological parameters.
The project aims to understand the mechanisms underlying cardiac arrhythmogenesis and defibrillation by running complex multi-scale computer models of the heart. Defibrillation by the timely application of strong electric shocks to the myocardium is the only effective therapy against cardiac arrhythmias. Due to the strong link existing between the upper limit of vulnerability (ULV) and the defibrillation threshold, a large body of research has aimed to unravel the mechanisms of cardiac vulnerability to electric shocks in order to better understand defibrillation failure. However, the contribution of transmural events to the mechanisms of shock-induced arrhythmogenesis is still unclear. This lack of understanding is in part due to the fact that current experimental techniques are unable to provide information regarding the transmural events that contribute to shock-induced arrhythmogenesis. Computer simulations using anatomically accurate ventricular models, in contrast, have the demonstrated ability to provide high-resolution insight into the electrical events that take place in the depth of the ventricular wall before, during and after the application of strong shocks.
Goal: Current computer models of the electrical activity within the heart treat the tissue as entirely homogeneous. However, there is experimental evidence stating that the walls of the ventricles can be divided into three distinct layers, containing different cell types (epicardial, midmyocardial and endocardial). The electrical differences between the three layers results in transmural dispersion in action potential duration (APD) which has been shown to contribute to arrhythmia induction in the heart. In order to make current computer models more realistic, the goal of this study is to use computer simulations to investigate how the inclusion of electrical heterogeneities affects cardiac vulnerability to electric shocks.
Methods: A 3D anatomically accurate, finite-element bidomain model of the rabbit ventricles was used. Transmural heterogeneities (TH) in ionic currents were incorporated based on experimental data. Computer simulations using the whole ventricular model were run on the National Grid Service (NGS), a distributed computing service that allows for large-scale computation. The ventricles were paced epicardially and after 7 paced beats, 8ms truncated exponential monophasic shocks of varying strength and coupling interval (CI) were applied via large planar external electrodes located at the boundaries of the perfusing bath. The vulnerable area (VA) was determined in the model with TH and compared to that obtained from a homogeneous APD model. The upper limit of vulnerability (ULV) is the highest shock strength above which no arrhythmia is induced and approximates the defibrillation threshold needed by clinicians.
Results: Simulations show that ventricles with heterogeneous APD profile are more vulnerable to electric shocks than those of homogeneous APD. This is due to both increased ULV, (from 26.7 to 30.5V/cm) and enlarged vulnerable area, (CIs from 110 to 170ms in the homogeneous case and from 120 to 190ms in case of TH). Analysis of results revealed that the elevation in ULV in the heterogeneous ventricles stems from increased transmural dispersion in postshock conduction velocity within the ventricular wall. In contrast, changes in VA are due to differences in transmural virtual electrode polarization (VEP) within the septum. For long CIs, 180-190ms, septal tissue in the homogeneous ventricles is strongly depolarized by the shock, and this results in propagation block shortly after shock end. In contrast, for the same CIs, the septum in the ventricles with TH exhibits a large excitable area, which facilitates the establishment of an intramural re-entrant circuit appearing as focal activity on the epicardium.
Conclusion: The use of computer models and simulations has allowed this project to investigate the effects of electrical heterogeneities on cardiac vulnerability to electric shocks in the 3D volume of the rabbit ventricles. Dispersion in postshock repolarisation within the LV wall plays a key role in the mechanisms underlying the increase in the ULV whereas shock-end polarisation within the septum determines the changes in the area of vulnerability.
We describe several two-sex population models exposed to a mild and long-lasting sexually transmitted disease. i.e. without disease-induced mortality and recovery. We modify these models to include non-reproductive groups and analyze their potential impact on the general population dynamics and of the disease in particular. The transmission of the disease is modeled through formation/separation of heterosexual couples assuming that one infected individual automatically infects his/her partner. We are interested in how the non-reproductive class may curb the growth of the infected group while keeping the healthy one at acceptable levels. A comparison with our previous results from one-sex models is also provided.
Tangle analysis has been applied successfully to study proteins which bind two segments of DNA and can knot and link circular DNA. We show how tangle analysis can be extended to model any stable protein-DNA complex.
Results: We have developed a computational algorithm to find the topological conformation of DNA bound within a protein complex. The algorithm uses an elementary invariant from knot theory called colorability to encode and search for possible DNA conformations. We apply this algorithm to analyze the experimental results of Pathania, Jayaram, and Harshey (Cell 2002). We show that the only topological DNA conformation bound by Mu transposase which is biologically likely is the five crossing solution found by Pathania {\it et al} (although other possibilities are discussed).
Conclusions: Our algorithm can be used to analyze the results of the experimental technique described in Pathania {\it et al} in order to determine the topological conformation of DNA bound within a stable protein-DNA complex.
In light of recent developments in rotavirus vaccine development and licensing, one contributing factor to assist policy makers on whether to add rotavirus vaccination to national immunization programs is to understand the costs and benefits associated with this type of policy decision. In Egypt, it is estimated that close to 1,813,550 will have become ill with rotavirus more and more than 3000 die due to rotavirus before reaching the age of five. To inform the decision makers at the Egyptian Ministry of Health and Population, a cost-benefit analysis, from the perspective of the MoHP, based on available local data from published and unpublished sources was conducted to evaluate the economic impact of introducing a rotavirus vaccine to the current national immunization schedule. A deterministic ODE-based model evaluating cross-immunity among wild-type strains and the vaccine strain is constructed to evaluate the vaccines impact on morbidity and mortality from rotavirus and to evaluate the effectiveness of different vaccination strategies.
Normally, the growth of neuronal processes is highly regulated, with temporally-controlled initial outgrowth followed by regular cycling between phases of extension and retraction as the processes navigate toward their targets. Ethanol delays the initial outgrowth of axons and dendrites, but the dynamics of subsequent growth are altered in rather specific ways [1]. We propose quantitative models that describe an axons random spatial displacement with and without a chemical gradient present. These models will be based on the observations that axons are lead by their growth cone towards a target in a random fashion. We believe this type of random response is caused by thermal fluctuations as well as the unpredictable intracellular dynamics of the growth cones response to spatial gradients. From assumptions based on experimental observations, we have assumed that Langevin dynamics is an adequate initial description, thus we propose a coupled SDE system. This system is studied in various ways, first via Ito calculus where I will derive statistical values from various coupled stochastic differential equation systems describing the neuron-axon-growth cone-gradient interaction. Other methods include the use of monte-carlo simulation techniques to obtain realizations of the biological process, as well as analyze the respective time evolution equation of the probability density function governing the SDE system. In comparing these math models to experimental data of a few cells we can see similar distributions in the lengthening behavior as well as the turning response (change in angle distribution). To further this study, we must employ multiple statistical time series analysis of extensive amounts of data for specific cells and compare these findings to our model simulations and solutions.
A recent randomized controlled trial shows a significant reduction in women-to-men transmission of HIV due to male circumcision. Such development calls for a rigorous mathematical study to ascertain the full impact of male circumcision in reducing HIV burden, especially in resource-poor nations where access to anti-retroviral drugs is limited. This poster presents a compartmental model for the transmission dynamics of HIV in a community where male circumcision is practiced. In addition to having a disease-free equilibrium, which is locally-asymptotically stable whenever a certain epidemiological threshold is less than unity, the model exhibits the phenomenon of backward bifurcation, where the disease-free equilibrium coexists with a stable endemic equilibrium when the threshold is less than unity. The implication of this result is that HIV may persist in the population even when the reproduction threshold is less than unity. Using partial data from South Africa, the study shows that male circumcision can prevent up to 150,000 cases and 9,400 deaths in the country within a year. Based on the estimate of circumcision efficacy from the randomized trial, this study shows that at least 60% circumcision coverage need to be attained to offer a realistic chance of effectively controlling HIV using circumcision alone. Further, it is shown that male circumcision will have positive impact (minimize HIV burden) if a certain quantity, defined as {it circumcision impact}, is positive; and will not otherwise.
Leukocyte rolling is a crucial part of the body's immune response and allows leukocytes from the bloodstream to reach sites of infection through receptor-ligand interactions with the vascular wall. Although many more advanced mathematical models exist to describe this type of cell adhesion, we consider a rigid cylinder in a shear flow near a flat wall as a two-dimensional toy model of a biological cell. The microscopic binding molecules are modelled as a continuous distribution of hookean springs and we investigate how their properties mediate macroscopic variables, such as the cell velocity. In contrast with other studies, each bond between the cell and the wall is allowed to tilt and form an angle with the normal to the wall. The cell, driven by the shear flow, is then no longer limited to a tank-treading motion but can also slide on the surface. Using asymptotic expansions we show that this extra degree of freedom in the model gives this mechanical system interesting properties like hys- teresis in the force-response relationship and a complete lift-off above a critical horizontal force. For a range of parameter values, we find a critical shear rate above which most of the bonds break and the cell moves freely over the wall. It will resume tank-treading only for a much lower shear rate. Albeit very simplified, the model sheds light on a friction mechanism that to our knowledge has not been studied in the context of cell adhesion and has macroscopic effects that could be observed experimentally.
Classic models of HIV dynamics describe infection using a single compartment with free mixing of virus and T cells. However, recent research has revealed significant differences in infection dynamics amongst the body's various spatially distinct lymphocyte compartments. To aid our understanding of the role the spatial structure of HIV infection, we build a differential-equation based model of T cell recirculation, homeostasis, and activation in the presence of HIV. This model helps us distinguish changes in T cell concentrations that are a direct consequence of virus-induced cell death from changes that are an indirect consequence recirculation. Model results are compared with experimental data, and the implications for HIV treatment are discussed.
At a genetic level humans and chimpanzees are closely related, with 98.7% of their DNA identical. It has long been speculated that many of the obvious differences between the two species are due to changes in regulatory sequences that control how genes are expressed (King and Wilson, 1975). A regulatory sequence is a short sequence of DNA (in vertebrates many are 6-9 nucleotides long) which is a binding site for transcription factors that promote or inhibit transcription of DNA to make proteins. Given what is known about transcription factor binding sites, this motivates the following probability question: Given a 1000 nucleotide region in our genome, how long does it take for a specified six to nine letter word to appear in that region in some individual? Stone and Wray (2001) computed 5,950 years as the answer for a given six letter word to appear in a 2000 nucleotide region. As MacArthur and Brookfield (2004) have already pointed out, there is a serious problem with this computation: they assumed that DNA sequences of different individuals are independent to get results for the whole population from their simulations for a single individual. In our paper, we show that if we do not assume independence, then for words of length 6, the average waiting time is 100,000 years. For words of length 8, the waiting time has mean 375,000 years if there is an almost match in the consensus sequence for the population, but if not the mean waiting time is 650 million years. Fortunately, in biological reality, the match to the target word does not have to be perfect for binding to occur. If we model this by saying that one mismatch is good enough, then the mean becomes about 60,000 years. We are currently extending this work to handle organisms with larger population sizes. The above results depend heavily on the fact that for the human effective population of size 10,000, double mutations between fixations are rare and we can ignore the possibility of triple mutations. This is false for organisms such as Drosophila (fruit flies), yeast, and E. coli. Here we consider the mathematical problem: How long does it take in a population of N (1,000,000 or more) diploid individuals until some individual has accumulated k specified mutations? Focusing on k = 2, 3, and 4, we do a similar analysis of the waiting time problem for these organisms.
A deterministic model for HIV, which accounts for its staged-progression property and viral load variability upon infection, is formulated and used to evaluate various treatment strategies, based on the use of anti-retroviral drugs. The strategies considered are to treat (i) every infected individual (universal strategy), (ii) individuals with low viral load only, (iii) individuals with high viral load (including those with AIDS symptoms) and (iv) people with clinical symptoms of AIDS only. The model, built via a refinement of a basic model, is analysed qualitatively, showing the existence of a globally-stable disease-free equilibrium whenever a certain threshold is less than unity and disease persistence whenever the threshold exceeds unity. The study shows that the universal strategy gives the highest reduction in the total number of cases, regardless of the treatment rate. Further, this strategy reduces more new cases than any of the others. This is followed by the high viral load strategy, then the AIDS-only strategy and finally the low viral load strategy. If treatment rates are low (perhaps due to limited supply of ARVs), the high viral load strategy, and the AIDS-only strategy avert more HIV-related mortality than any of the remaining strategies. For high treatment rates, the universal strategy averts more deaths than any of the other strategies. Finally, the study shows that treating individuals with low viral load could be detrimental to the community, since it may result in more deaths than the case where no treatment is used.
Three-dimensional collagen I gels are frequently used as a medium for studying cell motility, cancer invasion, and wound healing. The mechanical interactions between the cells and the gel play a critical role in directing cell motility, but these interactions are not well understood. Often, continuum models are used to describe gel mechanics, but these models typically ignore plastic deformations and fiber alignment. We present a discrete model of the gel, where each fiber is modeled individually and present comparisions between this model and experiment.
Many infections in both human and veterinary medicine remain difficult to manage clinically or prevent by vaccination. A successful immune response to infections must balance the pro-inflammatory effector actions with the protection of organ function. An anti-inflammatory component to the immune response is proven by the existence of host-derived immunosuppressive cytokines such as IL-10. The complexity of anti-inflammation regulation is revealed by the relatively recent description of CD4+ T regulatory cells. Although the importance of immunoregulation is appreciated in many experimental systems, the characteristics of infection that lead to changes in the balance of pro- vs. anti-inflammatory actions is largely unknown.
Working with immunologists, we have build up a sophisticated immune system in tissues involving the minimal set of key cells and signals known to play major role in immunity. Our PDE model not only successfully recapitulates major immune phenomena, but also test the well-known hypothesis in the immunology community and make sound predications in immune responses to various of infectious diseases.
In order to develop appropriate mathematical models of fluid-structure interaction which can be used in simulations of micro air vehicle (MAV) flight, it is crucial to understand how the elastic membrane component of a typical MAV wing will realistically behave under aerodynamic loads. Such wing structures are very thin and can exhibit large deflections, often of the same order of magnitude as the membrane cord length.
In this talk, we will outline a finite element procedure for computing the deflection of such membrane wings in response to aerodynamic pressure, where the structural model includes both geometric and material nonlinearity components.
In recent years, tumor-induced angiogenesis has become an important field of research because it represents a crucial step in the development of malignant tumors. The process is regulated by the interactions between various cell types such as endothelial cells (ECs) and macrophages, and by biochemical factors. These include angiogenic promoters such as vascular endothelial growth factor (VEGF) and inhibitors such as angiostatin.
In this poster we present a hybrid mathematical model in which cells are treated as discrete units in a continuum field of a chemoattractant(such as VEGF). The chemoattractant evolves according to a system of reaction-diffusion equations, whereas the discrete cells serve as sources/sinks in this continuum field. This discrete cell model for capillary tube formation of ECs incorporates a realistic model for signal transduction and VEGF production and release, and gives insights into the aggregation patterns and the factors that influence stream formation. In particular, it serves as a tool for investigating tumor-vessel signaling and the role of mechano-chemical interactions of the cells with the substratum. It therefore contributes towards the qualitative and quantitative understanding of angiogenesis, which may in turn lead to the development of new anti-angiogenic therapeutic strategies.
Bayesian network structure learning is a popular technique for attempting to model biological networks. This poster will highlight two related topics: 1) recent exploratory work which has investigated the appropriateness of different types of Bayesian network learning for modeling biological or biologically motivated datasets, and 2) preliminary work on relating the Bayesian network modeling approach to alternative mathematical modeling techniques.
The most damaging change during cancer progression is the growth of metastases. The protein RhoC GTPase was found to be crucial in that process in different cancers, particularly, in a highly aggressive form of breast cancer. RhoC is a molecular switch cycling between inactive (GDP-bound) and active (GTP-bound) states, tightly regulated by several regulatory proteins. We have developed a dual mathematical-experimental approach to understand this cycle and its deregulation in breast cancer cells in comparison with normal ones. A major impact of this work is to quantitatively predict the effects of drugs targeted against RhoC in cancer.
Several insulin therapies are widely in clinical use with the basic strategy that mimics insulin secretion in normal glucose-insulin endocrine metabolic regulatory system. In this paper, we model the glucose-insulin regulator system using systems of delay differential equations. We study the dynamics of the model both qualitatively and quantitatively. The analytical results show the existence and uniqueness of a stable periodic solution that corresponds to ultradian insulin secretion oscillations. Numerically we simulate the insulin administration based on our model. We confirm the effectiveness of the clinical applications of Replacement Therapy with basal-bolus insulin infusion for diabetes. The study will provide more efficient and effective algorithms for the treatments of diabetes mellitus in clinical applications.
This is a joint work with Drs. Jiaxu Li and Yang Kuang at the Department of Mathematics and Statistics at Arizona State University.
Many biological applications involve the dynamics of a large number of fibers immersed in a viscous biofluid. It is of significant importance to accurately predict the physical and biological behavior of the fiber-fluid interactions. We formulate a computational hydrodynamic model to attack such problems and make use of a numerical technique known as the immersed boundary method. We discuss the effects of fiber shapes, fluid viscosities and physical boundaries on the behavior of the fiber suspension and sedimentation.
Signaling through mitogen-activated protein kinase (MAPK) pathways is critical for cellular decisions to proliferate, differentiate, or undergo apoptosis. MAPK cascades are consist of several levels, where the activated kinase at each level phosphorylates the kinase at the next level down the cascade. Here we focus on one tier of it, a set of reactions of two-site MAPK phosphorylation and dephosphorylation. Using geometric singular perturbation theory, we prove that the system will always converge to some steady state, starting from almost every initial condition. The same result holds for other time scale separated systems, provided the limiting systems are strongly monotone.
The transcription factor p53 plays a central role in maintaining genomic integrity. Recent experiments have shown that p53 protein level rises and falls in distinct pulses in response to DNA damage. The amplitudes of and intervals between pulses seem to be independent of the extent of damage, but extensive or irreparable damage may lead to sustained pulsing of p53. Identifying the molecular mechanisms responsible for such interesting behavior is an important and challenging problem. Our work describes four dual-feedback mechanisms that combine both positive and negative feedback loops, which have been identified in the signaling network responsible for p53 regulation. Mathematical models of all four mechanisms are analyzed to determine if they are consistent with experimental observations and to characterize subtle differences among the possible mechanisms. In addition, a novel molecular mechanism is proposed whereby p53 pulses may induce, at first, cell cycle arrest and, if sustained, cell death. The proposal accounts for basic features of p53-mediated responses to DNA damage and suggests new experiments to probe the dynamics of p53 signaling.
