Novel recurrent neural network for modelling biological networks: Oscillatory p53 interaction dynamics

Understanding the control of cellular networks consisting of gene and protein interactions and their emergent properties is a central activity of Systems Biology research. For this, continuous, discrete, hybrid, and stochastic methods have been proposed. Currently, the most common approach to modelling accurate temporal dynamics of networks is ordinary differential equations (ODE). However, critical limitations of ODE models are difficulty in kinetic parameter estimation and numerical solution of a large number of equations, making them more suited to smaller systems. In this article, we introduce a novel recurrent artificial neural network (RNN) that addresses above limitations and produces a continuous model that easily estimates parameters from data, can handle a large number of molecular interactions and quantifies temporal dynamics and emergent systems properties. This RNN is based on a system of ODEs representing molecular interactions in a signalling network. Each neuron represents concentration change of one molecule represented by an ODE. Weights of the RNN correspond to kinetic parameters in the system and can be adjusted incrementally during network training. The method is applied to the p53-Mdm2 oscillation system – a crucial component of the DNA damage response pathways activated by a damage signal. Simulation results indicate that the proposed RNN can successfully represent the behaviour of the p53-Mdm2 oscillation system and solve the parameter estimation problem with high accuracy. Furthermore, we presented a modified form of the RNN that estimates parameters and captures systems dynamics from sparse data collected over relatively large time steps. We also investigate the robustness of the p53-Mdm2 system using the trained RNN under various levels of parameter perturbation to gain a greater understanding of the control of the p53-Mdm2 system. Its outcomes on robustness are consistent with the current biological knowledge of this system. As more quantitative data become available on individual proteins, the RNN would be able to refine parameter estimation and mapping of temporal dynamics of individual signalling molecules as well as signalling networks as a system. Moreover, RNN can be used to modularise large signalling networks.     

Estimating a predator-prey dynamical model with the parameter cascades method

Summary .   Ordinary differential equations (ODEs) are widely used in ecology to describe the dynamical behavior of systems of interacting populations. However, systems of ODEs rarely provide quantitative solutions that are close to real field observations or experimental data because natural systems are subject to environmental and demographic noise and ecologists are often uncertain about the correct parameterization. In this article we introduce "parameter cascades" as an improved method to estimate ODE parameters such that the corresponding ODE solutions fit the real data well. This method is based on the modified penalized smoothing with the penalty defined by ODEs and a generalization of profiled estimation, which leads to fast estimation and good precision for ODE parameters from noisy data. This method is applied to a set of ODEs originally developed to describe an experimental predator–prey system that undergoes oscillatory dynamics. The new parameterization considerably improves the fit of the ODE model to the experimental data sets. At the same time, our method reveals that important structural assumptions that underlie the original ODE model are essentially correct. The mathematical formulations of the two nonlinear interaction terms (functional responses) that link the ODEs in the predator–prey model are validated by estimating the functional responses nonparametrically from the real data. We suggest two major applications of "parameter cascades" to ecological modeling: It can be used to estimate parameters when original data are noisy, missing, or when no reliable priori estimates are available; it can help to validate the structural soundness of the mathematical modeling approach.     

Metal reduction kinetics in Shewanella

MOTIVATION: Metal reduction kinetics have been studied in cultures of dissimilatory metal reducing bacteria which include the Shewanella oneidensis strain MR-1. Estimation of system parameters from time-series data faces obstructions in the implementation depending on the choice of the mathematical model that captures the observed dynamics. The modeling of metal reduction is often based on Michaelis-Menten equations. These models are often developed using initial in vitro reaction rates and seldom match with in vivo reduction profiles. RESULTS: For metal reduction studies, we propose a model that is based on the power law representation that is effectively applied to the kinetics of metal reduction. The method yields reasonable parameter estimates and is illustrated with the analysis of time-series data that describes the dynamics of metal reduction in S.oneidensis strain MR-1. In addition, mixed metal studies involving the reduction of Uranyl (U(VI)) to the relatively insoluble tetravalent form (U(IV)) by S. alga strain (BR-Y) were studied in the presence of environmentally relevant iron hydrous oxides. For mixed metals, parameter estimation and curve fitting are accomplished with a generalized least squares formulation that handles systems of ordinary differential equations and is implemented in Matlab. It consists of an optimization algorithm (Levenberg-Marquardt, LSQCURVEFIT) and a numerical ODE solver. Simulation with the estimated parameters indicates that the model captures the experimental data quite well. The model uses the estimated parameters to predict the reduction rates of metals and mixed metals at varying concentrations. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.     

Understanding the transmission dynamics of respiratory syncytial virus using multiple time series and nested models

The nature and role of re-infection and partial immunity are likely to be important determinants of the transmission dynamics of human respiratory syncytial virus (hRSV). We propose a single model structure that captures four possible host responses to infection and subsequent reinfection: partial susceptibility, altered infection duration, reduced infectiousness and temporary immunity (which might be partial). The magnitude of these responses is determined by four homotopy parameters, and by setting some of these parameters to extreme values we generate a set of eight nested, deterministic transmission models. In order to investigate hRSV transmission dynamics, we applied these models to incidence data from eight international locations. Seasonality is included as cyclic variation in transmission. Parameters associated with the natural history of the infection were assumed to be independent of geographic location, while others, such as those associated with seasonality, were assumed location specific. Models incorporating either of the two extreme assumptions for immunity (none or solid and lifelong) were unable to reproduce the observed dynamics. Model fits with either waning or partial immunity to disease or both were visually comparable. The best fitting structure was a lifelong partial immunity to both disease and infection. Observed patterns were reproduced by stochastic simulations using the parameter values estimated from the deterministic models.     

Spatial models of virus-immune dynamics

To date, the majority of theoretical models describing the dynamics of infectious diseases in vivo are based on the assumption of well-mixed virus and cell populations. Because many infections take place in solid tissues, spatially structured models represent an important step forward in understanding what happens when the assumption of well-mixed populations is relaxed. Here, we explore models of virus and virus-immune dynamics where dispersal of virus and immune effector cells was constrained to occur locally. The stability properties of our spatial virus-immune dynamics models remained robust under almost all biologically plausible dispersal schemes, regardless of their complexity. The various spatial dynamics were compared to the basic non-spatial dynamics and important differences were identified: When space was assumed to be homogeneous, the dynamics generated by non-spatial and spatially structured models differed substantially at the peak of the infection. Thus, non-spatial models may lead to systematic errors in the estimates of parameters underlying acute infection dynamics. When space was assumed to be heterogeneous, spatial coupling not only changed the equilibrium properties of the uncoupled populations but also equalized the dynamics and thereby reduced the likelihood of dynamic elimination of the infection. In line with experimental and clinical observations, long-lasting oscillation periods were virtually absent. When source-sink dynamics were considered, the long-term outcome of the infection depended critically on the degree of spatial coupling. The infection collapsed when emigration from source sites became too large. Finally, we discuss the implications of spatially structured models on medical treatment of infectious diseases, and note that a huge gap exists in data accurately describing infection dynamics in solid tissues.     

HIV model parameter estimates from interruption trial data including drug efficacy and reservoir dynamics

Mathematical models based on ordinary differential equations (ODE) have had significant impact on understanding HIV disease dynamics and optimizing patient treatment. A model that characterizes the essential disease dynamics can be used for prediction only if the model parameters are identifiable from clinical data. Most previous parameter identification studies for HIV have used sparsely sampled data from the decay phase following the introduction of therapy. In this paper, model parameters are identified from frequently sampled viral-load data taken from ten patients enrolled in the previously published AutoVac HAART interruption study, providing between 69 and 114 viral load measurements from 3-5 phases of viral decay and rebound for each patient. This dataset is considerably larger than those used in previously published parameter estimation studies. Furthermore, the measurements come from two separate experimental conditions, which allows for the direct estimation of drug efficacy and reservoir contribution rates, two parameters that cannot be identified from decay-phase data alone. A Markov-Chain Monte-Carlo method is used to estimate the model parameter values, with initial estimates obtained using nonlinear least-squares methods. The posterior distributions of the parameter estimates are reported and compared for all patients.     

Numerical discretization-based estimation methods for ordinary differential equation models via penalized spline smoothing with applications in biomedical research

Differential equations are extensively used for modeling dynamics of physical processes in many scientific fields such as engineering, physics, and biomedical sciences. Parameter estimation of differential equation models is a challenging problem because of high computational cost and high-dimensional parameter space. In this article, we propose a novel class of methods for estimating parameters in ordinary differential equation (ODE) models, which is motivated by HIV dynamics modeling. The new methods exploit the form of numerical discretization algorithms for an ODE solver to formulate estimating equations. First, a penalized-spline approach is employed to estimate the state variables and the estimated state variables are then plugged in a discretization formula of an ODE solver to obtain the ODE parameter estimates via a regression approach. We consider three different order of discretization methods, Euler's method, trapezoidal rule, and Runge-Kutta method. A higher-order numerical algorithm reduces numerical error in the approximation of the derivative, which produces a more accurate estimate, but its computational cost is higher. To balance the computational cost and estimation accuracy, we demonstrate, via simulation studies, that the trapezoidal discretization-based estimate is the best and is recommended for practical use. The asymptotic properties for the proposed numerical discretization-based estimators are established. Comparisons between the proposed methods and existing methods show a clear benefit of the proposed methods in regards to the trade-off between computational cost and estimation accuracy. We apply the proposed methods t an HIV study to further illustrate the usefulness of the proposed approaches.     

The minimal model of the hypothalamic–pituitary–adrenal axis

This paper concerns ODE modeling of the hypothalamic–pituitary– adrenal axis (HPA axis) using an analytical and numerical approach, combined with biological knowledge regarding physiological mechanisms and parameters. The three hormones, CRH, ACTH, and cortisol, which interact in the HPA axis are modeled as a system of three coupled, nonlinear differential equations. Experimental data shows the circadian as well as the ultradian rhythm. This paper focuses on the ultradian rhythm. The ultradian rhythm can mathematically be explained by oscillating solutions. Oscillating solutions to an ODE emerges from an unstable fixed point with complex eigenvalues with a positive real parts and a non-zero imaginary parts. The first part of the paper describes the general considerations to be obeyed for a mathematical model of the HPA axis. In this paper we only include the most widely accepted mechanisms that influence the dynamics of the HPA axis, i.e. a negative feedback from cortisol on CRH and ACTH. Therefore we term our model the minimal model. The minimal model, encompasses a wide class of different realizations, obeying only a few physiologically reasonable demands. The results include the existence of a trapping region guaranteeing that concentrations do not become negative or tend to infinity. Furthermore, this treatment guarantees the existence of a unique fixed point. A change in local stability of the fixed point, from stable to unstable, implies a Hopf bifurcation; thereby, oscillating solutions may emerge from the model. Sufficient criteria for local stability of the fixed point, and an easily applicable sufficient criteria guaranteeing global stability of the fixed point, is formulated. If the latter is fulfilled, ultradian rhythm is an impossible outcome of the minimal model and all realizations thereof. The second part of the paper concerns a specific realization of the minimal model in which feedback functions are built explicitly using receptor dynamics. Using physiologically reasonable parameter values, along with the results of the general case, it is demonstrated that un-physiological values of the parameters are needed in order to achieve local instability of the fixed point. Small changes in physiologically relevant parameters cause the system to be globally stable using the analytical criteria. All simulations show a globally stable fixed point, ruling out periodic solutions even when an investigation of the ‘worst case parameters’ is performed.     

Estimation of the Distribution of Infection Times Using Longitudinal Serological Markers of HIV: Implications for the Estimation of HIV Incidence

Summary In the last decade, interest has been focused on human immunodeficiency virus (HIV) antibody assays and testing strategies that could distinguish recent infections from established infection in a single serum sample. Incidence estimates are obtained by using the relationship between prevalence, incidence, and duration of recent infection (window period). However, recent works demonstrated limitations of this approach due to the use of an estimated mean “window period.” We propose an alternative approach that consists in estimating the distribution of infection times based on serological marker values at the moment when the infection is first discovered. We propose a model based on the repeated measurements of virological markers of seroconversion for the marker trajectory. The parameters of the model are estimated using data from a cohort of HIV‐infected patients enrolled during primary infection. This model can be used for estimating the distribution of infection times for newly HIV diagnosed subjects reported in a HIV surveillance system. An approach is proposed for estimating HIV incidence from these results.     

A dynamical model of human immune response to influenza A virus infection

We present a simplified dynamical model of immune response to uncomplicated influenza A virus (IAV) infection, which focuses on the control of the infection by the innate and adaptive immunity. Innate immunity is represented by interferon-induced resistance to infection of respiratory epithelial cells and by removal of infected cells by effector cells (cytotoxic T-cells and natural killer cells). Adaptive immunity is represented by virus-specific antibodies. Similar in spirit to the recent model of Bocharov and Romanyukha [1994. Mathematical model of antiviral immune response. III. Influenza A virus infection. J. Theor. Biol. 167, 323-360], the model is constructed as a system of 10 ordinary differential equations with 27 parameters characterizing the rates of various processes contributing to the course of disease. The parameters are derived from published experimental data or estimated so as to reproduce available data about the time course of IAV infection in a na?ve host. We explore the effect of initial viral load on the severity and duration of the disease, construct a phase diagram that sheds insight into the dynamics of the disease, and perform sensitivity analysis on the model parameters to explore which ones influence the most the onset, duration and severity of infection. To account for the variability and speed of adaptation of the adaptive response to a particular virus strain, we introduce a variable that quantifies the antigenic compatibility between the virus and the antibodies currently produced by the organism. We find that for small initial viral load the disease progresses through an asymptomatic course, for intermediate value it takes a typical course with constant duration and severity of infection but variable onset, and for large initial viral load the disease becomes severe. This behavior is robust to a wide range of parameter values. The absence of antibody response leads to recurrence of disease and appearance of a chronic state with nontrivial constant viral load.     

Estimating kinetic parameters from HIV primary infection data through the eyes of three different mathematical models

The dynamics of HIV-1 infection consist of three distinct phases starting with primary infection, then latency and finally AIDS or drug therapy. In this paper we model the dynamics of primary infection and the beginning of latency. We show that allowing for time delays in the model better predicts viral load data when compared to models with no time delays. We also find that our model of primary infection predicts the turnover rates for productively infected T cells and viral totals to be much longer than compared to data from patients receiving anti-viral drug therapy. Hence the dynamics of the infection can change dramatically from one stage to the next. However, we also show that with the data available the results are highly sensitive to the chosen model. We compare the results using analysis and Monte Carlo techniques for three different models and show how each predicts rather dramatic differences between the fitted parameters. We show, using a chi(2) test, that these differences between models are statistically significant and using a jackknifing method, we find the confidence intervals for the parameters. These differences in parameter estimations lead to widely varying conclusions about HIV pathogenesis. For instance, we find in our model with time delays the existence of a Hopf bifurcation that leads to sustained oscillations and that these oscillations could simulate the rapid turnover between viral strains and the appropriate CTL response necessary to control the virus, similar to that of a predator-prey type system.     

Lessons on pattern formation from planet WATOR

It is well known that if reacting species experience unequal diffusion rates, then dynamics that lead to a constant steady state in a "well-mixed" environment can in a spatial setting lead to interesting patterns. In this paper, we focus on complementary pattern formation mechanisms that operate even when the diffusion rates are equal. In particular, we can say that when the mean-field ODE has an attracting periodic orbit then the stochastic spatial model will have large-scale spatial structures in equilibrium. We explore this mechanism in depth through the dynamics of the simulator WATOR.     

Robust estimation for ordinary differential equation models

Applied scientists often like to use ordinary differential equations (ODEs) to model complex dynamic processes that arise in biology, engineering, medicine, and many other areas. It is interesting but challenging to estimate ODE parameters from noisy data, especially when the data have some outliers. We propose a robust method to address this problem. The dynamic process is represented with a nonparametric function, which is a linear combination of basis functions. The nonparametric function is estimated by a robust penalized smoothing method. The penalty term is defined with the parametric ODE model, which controls the roughness of the nonparametric function and maintains the fidelity of the nonparametric function to the ODE model. The basis coefficients and ODE parameters are estimated in two nested levels of optimization. The coefficient estimates are treated as an implicit function of ODE parameters, which enables one to derive the analytic gradients for optimization using the implicit function theorem. Simulation studies show that the robust method gives satisfactory estimates for the ODE parameters from noisy data with outliers. The robust method is demonstrated by estimating a predator-prey ODE model from real ecological data.     

Modeling and Estimation of Kinetic Parameters and Replicative Fitness of HIV-1 from Flow-Cytometry-Based Growth Competition Experiments

Growth competition assays have been developed to quantify the relative fitness of HIV-1 mutants. In this article, we develop mathematical models to describe viral/cellular dynamic interactions in the assay system from which the competitive fitness indices or parameters are defined. In our previous HIV-viral fitness experiments, the concentration of uninfected target cells was assumed to be constant (Wu et al. 2006). But this may not be true in some experiments. In addition, dual infection may frequently occur in viral fitness experiments and may not be ignorable. Here, we relax these two assumptions and extend our earlier viral fitness model (Wu et al. 2006). The resulting models then become nonlinear ODE systems for which closed-form solutions are not achievable. In the new model, the viral relative fitness is a function of time since it depends on the target cell concentration. First, we studied the structure identifiability of the nonlinear ODE models. The identifiability analysis showed that all parameters in the proposed models are identifiable from the flow-cytometry-based experimental data that we collected. We then employed a global optimization approach (the differential evolution algorithm) to directly estimate the kinetic parameters as well as the relative fitness index in the nonlinear ODE models using nonlinear least square regression based on the experimental data. Practical identifiability was investigated via Monte Carlo simulations.     

Uncertainties regarding dengue modeling in Rio de Janeiro, Brazil

Dengue fever is currently the most important arthropod-borne viral disease in Brazil. Mathematical modeling of disease dynamics is a very useful tool for the evaluation of control measures. To be used in decision-making, however, a mathematical model must be carefully parameterized and validated with epidemiological and entomological data. In this work, we developed a simple dengue model to answer three questions: (i) which parameters are worth pursuing in the field in order to develop a dengue transmission model for Brazilian cities; (ii) how vector density spatial heterogeneity influences control efforts; (iii) with a degree of uncertainty, what is the invasion potential of dengue virus type 4 (DEN-4) in Rio de Janeiro city. Our model consists of an expression for the basic reproductive number (R0) that incorporates vector density spatial heterogeneity. To deal with the uncertainty regarding parameter values, we parameterized the model using a priori probability density functions covering a range of plausible values for each parameter. Using the Latin Hypercube Sampling procedure, values for the parameters were generated. We conclude that, even in the presence of vector spatial heterogeneity, the two most important entomological parameters to be estimated in the field are the mortality rate and the extrinsic incubation period. The spatial heterogeneity of the vector population increases the risk of epidemics and makes the control strategies more complex. At last, we conclude that Rio de Janeiro is at risk of a DEN-4 invasion. Finally, we stress the point that epidemiologists, mathematicians, and entomologists need to interact more to find better approaches to the measuring and interpretation of the transmission dynamics of arthropod-borne diseases.     

Steady state analysis of influx and transbilayer distribution of ergosterol in the yeast plasma membrane

Mathematical modeling is now frequently used in outbreak investigations to understand underlying mechanisms of infectious disease dynamics, assess patterns in epidemiological data, and forecast the trajectory of epidemics. However, the successful application of mathematical models to guide public health interventions lies in the ability to reliably estimate model parameters and their corresponding uncertainty. Here, we present and illustrate a simple computational method for assessing parameter identifiability in compartmental epidemic models. We describe a parametric bootstrap approach to generate simulated data from dynamical systems to quantify parameter uncertainty and identifiability. We calculate confidence intervals and mean squared error of estimated parameter distributions to assess parameter identifiability. To demonstrate this approach, we begin with a low-complexity SEIR model and work through examples of increasingly more complex compartmental models that correspond with applications to pandemic influenza, Ebola, and Zika. Overall, parameter identifiability issues are more likely to arise with more complex models (based on number of equations/states and parameters). As the number of parameters being jointly estimated increases, the uncertainty surrounding estimated parameters tends to increase, on average, as well. We found that, in most cases, R0 is often robust to parameter identifiability issues affecting individual parameters in the model. Despite large confidence intervals and higher mean squared error of other individual model parameters, R0 can still be estimated with precision and accuracy. Because public health policies can be influenced by results of mathematical modeling studies, it is important to conduct parameter identifiability analyses prior to fitting the models to available data and to report parameter estimates with quantified uncertainty. The method described is helpful in these regards and enhances the essential toolkit for conducting model-based inferences using compartmental dynamic models.     

The Dynamics of a SEIR–SIRC Antigenic Drift Influenza Model

We consider the dynamics of an influenza model with antigenic drift mechanism. Antigenic drift is an antigen mutation on the skin surface of the influenza virus that do not produce a new virus strain. The mutation produces the same virus but with slightly different antigens that cannot be recognized by the immune receptors formed by the previous infection. There are some type of influenza that involve the interaction between two populations such as human and animal. In this paper, we construct an influenza model with antigenic drift mechanism on the human population that has interaction with the animal population. The animal population is assumed to follow the SEIR epidemic model. Our model is motivated by the fact that some of the influenza cases in human come from the animal such as the swine and the avian. The transmission parameter that shows number of contact between the susceptible human and the infectious animals are important to study. The parameter plays an important role to detect the cycle of infection of the disease. The other important parameters are the seasonality degree, which shows the pathogen appearance and disappearance via annual migration, and the infection rate on the human population. We employ the bifurcation theory to analyze the behavior of the system and to detect the cycle of infection types when the parameters values are varied.     

Effect of diffusion and spatially varying predation risk on the dynamics and equilibrium density of a predator-prey system

Starting from natural planktonic systems, we present a new mechanism involving spatial heterogeneity, and develop a new spatial structure model of planktonic predation systems. Firstly, the effect of diffusion on the dynamics of the system is investigated. We find that diffusion of only prey or both prey and predator between different patches with different predation risk may stabilize the dynamics, depending on the flow rate. Only a medium flow rate can lead to the stability of the system. Too large a rate can cause the system to approach the non-spatial limit case of a well-mixed system. Too large a rate can cause the system to approach the non-spatial limit case as a well-mixed system, which is characterized by its strongly oscillatory dynamics. When only prey diffuse, the smaller the parameter f (the proportion of the patchy volume with larger predation risk to the total volume), the more stable the system. If both populations can diffuse, however, only medium and very small f values may stabilize the system. Also, the response of the spatially averaged equilibrium densities of the system to the increasing of the flow rate is examined. With increasing flow rate, the spatial-averaged equilibrium density of prey decreases, while that of predator depends on which species can diffuse. For the case of prey diffusion only, it first remains unchanged and then slightly decreases, while it increases for the case of combinations as the flow rate increases. Our results are, qualitatively, determined by the spatially heterogeneous mechanism that we propose, and further regulated by top-down forces. Of practical importance, the results reported here indicate that which species diffuse plays a key role in the ways in which diffusion influences the dynamics and the spatial-average equilibrium densities of the system responses to the flow rate's increasing.