A mean-field approach for modeling the propagation of perturbations in biochemical reaction networks

Abstract

Often, the time evolution of a biochemical reaction network is crucial for determining the effects of combining multiple pharmaceuticals. Here we illustrate a mathematical framework for modeling the dominant temporal behaviour of a complicated molecular pathway or biochemical reaction network in response to an arbitrary perturbation, such as resulting from the administration of a therapeutic agent. The method enables the determination of the temporal evolution of a target protein as the perturbation propagates through its regulatory network. The mathematical approach is particularly useful when the experimental data that is available for characterizing or parameterizing the regulatory network is limited or incomplete. To illustrate the method, we consider the examples of the regulatory networks for the target proteins c-Myc and Chop, which play an important role in venetoclax resistance in acute myeloid leukemia. First we show how the networks that regulate each target protein can be reduced to a mean-field model by identifying the distinct effects that groups of proteins in the regulatory network have on the target protein. Then we show how limited protein-level data can be used to further simplify the mean-field model to pinpoint the dominant effects of the network perturbation on the target protein. This enables a further reduction in the number of parameters in the model. The result is an ordinary differential equation model that captures the temporal evolution of the expression of a target protein when one or more proteins in its regulatory network have been perturbed. Finally, we show how the dominant effects predicted by the mathematical model agree with RNA sequencing data for the regulatory proteins comprising the molecular network, despite the model not having a priori knowledge of this data. Thus, while the approach gives a simplified model for the expression of the target protein, it allows for the interpretation of the effects of the perturbation on the regulatory network itself. This method can be easily extended to sets of target proteins to model components of a larger systems biology model, and provides an approach for partially integrating RNA sequencing data and protein expression data. Moreover, it is a general approach that can be used to study drug effects on specific protein(s) in any disease or condition.

Introduction

Systems biology has emerged in recent years as a powerful method for quantitative modeling of biological phenomena. The approach has gained attention because it is generic and it enables the integration of experimental data with mathematical and computational modeling to understand and make predictions about complex biological processes. Using a systems biology approach, one can model virtually any type of interacting system, including protein and gene networks, metabolic networks, cellular interactions, and tissue and organism-level interactions. The effectiveness of the method has contributed to its popularity for addressing questions in a number of different biological settings, including in molecular pathways underlying cancers and other diseases Barillot et al. (2012); Chuang et al. (2010); Ebhardt et al. (2015); Faratian et al. (2009); Greenblum et al. (2012); Hornberg et al. (2006); Kreeger and Lauffenburger (2010); Loscalzo and Barabasi (2011); Mani et al. (2008); Nicholson (2006); Parikshak et al. (2015); Perou and Børresen-Dale (2011); Yarden and Pines (2012), in drug discovery and target identification Arrell and Terzic (2010); Azmi et al. (2010); Berg et al. (2005); Bugrim et al. (2004); Butcher (2005); Butcher et al. (2004); Keskin et al. (2007); Zhu et al. (2008), and in overcoming drug resistance Galluzzi et al. (2014); Wang et al. (2013).

However, one limitation that arises in systems biology approaches, particularly in the context of molecular pathways and biochemical reaction networks, such as gene and protein regulatory networks, is that they often include several kinetic parameters that cannot be determined independently of experimental time series data. The limited availability of quantitative experimental data then leads to a trade-off in the development of the model. On one hand, the goal of a systems biology approach is to develop a mathematical model that is complex enough to capture the observed biological phenomena and to answer questions related to the corresponding protein or gene network. On the other hand, increasing the complexity of the model without sufficient experimental data leads to a model with poor predictive power Ay and Arnosti (2011). In practice, fitting the outputs of the systems biology model to a limited experimental data set can lead to non-unique parameter sets or subsets of parameters that give fits of similar quality. This results in difficulties in interpreting the meanings of the parameters in the underlying molecular pathway or biochemical reaction network.

There are other methods for modeling the temporal evolution of gene regulatory networks, such as boolean models, that require only qualitative experimental data. However, while these types of models are easily extendable to large-scale systems, they tend to produce only qualitative results Ay and Arnosti (2011). On the other hand, several quantitative model reduction techniques have been developed to overcome the challenges posed by large-scale systems biology models. Such methods were recently reviewed in Snowden et al. (2017) and can be roughly categorized into timescale exploitation approaches Briggs and Haldane (1925); Choi et al. (2008); Debussche and Temam (1991); Härdin et al. (2009); Klonowski (1983); Kokotović (1984); Lam (1985); Maas and Pope (1992); Michaelis and Menten (1913); Noel et al. (2012); Prescott and Papachristodoulou (2014); Schneider and Wilhelm (2000); Surovtsova et al. (2009); Tikhonov (1952); Vejchodskỳ et al. (2014); West et al. (2015); Zobeley et al. (2005), reduction via sensitivity analysis Apri et al. (2012); Maurya et al. (2005a); Turanyi et al. (1989); Zhang and Goutsias (2010); Zi (2011), optimisation methods Anderson et al. (2011); Danø et al. (2006); Hangos et al. (2013); Maurya, Bornheimer, Venkatasubramanian, Subramaniam, 2009, Maurya, Scott, Venkatasubramanian, Subramaniam, 2005; Prescott and Papachristodoulou (2012); Taylor et al. (2008), lumping Dokoumetzidis and Aarons (2009); Koschorreck et al. (2007); Kuo and Wei (1969); Li and Rabitz (1990); Li, Rabitz, Tóth, 1994, Li, Tomlin, Rabitz, Tóth, 1994; Sunnåker, Cedersund, Jirstrand, 2011, Sunnåker, Schmidt, Jirstrand, Cedersund, 2010; Tomlin et al. (1994); Wei and Kuo (1969), and singular value decomposition-based approaches Hahn and Edgar (2000); Hardin and van Schuppen (2006); Lall et al. (2002); Liebermeister et al. (2005); Moore (1981); Sootla and Anderson (2014). Other approaches that do not strictly fit into these categories have also been developed for biochemical reaction networks, for example, based on iterative methods Apri et al. (2014); Maiwald et al. (2016); Quaiser et al. (2011); Rao et al. (2013) and differential geometry Transtrum, Qiu, 2014, Transtrum, Qiu, 2016. Each approach has its advantages, limitations, and distinct ranges of applicability Snowden et al. (2017). The optimal technique for a particular system depends on the context of the problem at hand, including the nature of the full-scale model and available experimental data and the goals of the modeller.

Model reduction approaches have also been developed that draw inspiration from mean-field theory in physics and probability theory Chaikin et al. (1995); Curie (1895); Weiss (1907), in which high-dimensional stochastic models are reduced to a simpler approximate model by averaging over the degrees of freedom in the system. This enables the reduction of a many-body problem to a one-body problem by approximating individual interactions by a single averaged effect. In the context of gene regulatory networks, mean-field approaches have been developed to construct simplified models of gene transcription Andrecut and Kauffman (2006); Cao and Grima (2018); Sturrock et al. (2015). However, in the case of  Andrecut and Kauffman (2006) and Sturrock et al. (2015), the resulting models take the form of systems of coupled ordinary differential equations with several kinetic parameters and are not particularly suitable for modeling large-scale networks. Indeed, the existing mean-field approaches were applied to the study of small-scale networks with a limited number of feedback loops and are well-suited for modeling gene regulatory networks that only include protein-promoter type interactions.

While several approaches have been developed to construct simplified models of gene regulatory networks, another common issue that arises in the mathematical modeling of genetic pathways is that, in many cases, transcript levels alone are not sufficient to predict protein levels Li and Biggin (2015); Liu et al. (2016). This is because RNA expression often does not correlate with protein expression for a number of reasons, including the availability of amino acids, the formation of higher order protein complexes, and post-translational modifications, such as phosphorylation and proteolytic events. Despite these biological constraints, often it is dubiously assumed, for simplicity, that the rate of protein translation is directly proportional to RNA concentration. Indeed, this is the case for the mean-field approaches discussed above.

In this work, we introduce a simplistic mean-field approach to model the temporal evolution of protein expression levels which partially overcomes the limitations and challenges faced by other modeling techniques. The method relies on the decomposition of a complicated molecular pathway into a simplified uncoupled model of the net effects exerted by the interaction network on a target protein. The method can be easily scaled to large interaction networks and can be generalized to encompass multiple target proteins and arbitrary types of perturbations to the network. Intriguingly, the dominant effects identified by the mean-field model, using only the target protein expression levels, are supported by RNA sequencing data of the regulatory proteins in the pathway, despite the model having no a priori knowledge of this data. Thus this approach offers interpretability to the modeling results and provides a stepping stone toward the integration of protein level data with gene expression data.

The manuscript is organized as follows. In Section 2, we discuss the relevant biological background and protein interactions for the examples considered in this work. In Section 3, we describe the experimental methods and the modeling approach, then present the results and discussion in Section 4. Finally, in Section 5 we give some concluding remarks.