Abstract
This protocol illustrates a pipeline for modeling the nonlinear behavior of intracellular signaling pathways. At fixed spatial points, nonlinear signaling dynamics are described by ordinary differential equations (ODEs). At constant parameters, these ODEs may have multiple attractors, such as multiple steady states or limit cycles. Standard optimization procedures finetune the parameters for the system trajectories localized within the basin of attraction of only one attractor, usually a stable steady state. The suggested protocol samples the parameter space and captures the overall dynamic behavior by analyzing the number and stability of steady states and the shapes of the assembly of nullclines, which are determined as projections of quasisteadystate trajectories into different 2D spaces of system variables. Our pipeline allows identifying main qualitative features of the model behavior, perform bifurcation analysis, and determine the borders separating the different dynamical regimes within the assembly of 2D parametric planes. Partial differential equation (PDE) systems describing the nonlinear spatiotemporal behavior are derived by coupling fixed point dynamics with species diffusion.
Keywords: Cell signaling, Nonlinear dynamics, Multistability, Oscillations, Bifurcations, Ordinary and partial differential equations
Background
Here, we present a protocol for computational analysis of the nonlinear dynamic behavior of signaling networks and their transitions between different dynamic regimes. The dynamics of spatially homogenous biochemical systems are described by ordinary differential equations (ODE), whereas the spatiotemporal dynamics are described by partial differential equation (PDE) systems. Knowledge about ODE and PDE is needed as a prerequisite for understanding this protocol. Our computational protocol analyzes the systems’ dynamics at fixed spatial points, and it considers the nonlinear spatiotemporal behavior by coupling fixed point dynamics with species diffusion (Tsyganov et al., 2012; BoladoCarrancio et al., 2020). A variable x_{i} that is the concentration or activity of each network node i depends on the time (t) in an ODE system and on both the time and the spatial coordinates in a PDE system. The described protocol operates with the already established interaction topology of the signaling network under study; that is, for each node i it is known what nodes activate and/or inhibit it. Mathematically, the network topology is determined by the signed incidence matrix of the ODE system,
If node j activates or inhibits node i (∂f_{i}/∂x_{j}>0 or ∂f_{i}/∂x_{j},<0, respectively), the element (i,j)of the signed incidence matrix equals 1 or 1, respectively, and it equals zero if node j does not affect node i.
Traditional modeling pipelines for ODE systems use parameter optimization procedures, which rely on the fitting of simulated trajectories to experimental data points, and then study the behavior of a calibrated model (Glover et al., 2000; Maiwald and Timmer, 2008; Kreutz and Timmer, 2009; Penas et al., 2015; Thomas et al., 2016; Degasperi et al., 2017; Mitra et al., 2018 and 2019). These pipelines produce satisfactory results when the system under study has a single stable steady state or when the system trajectories are within the basin of attraction of a stable steady state for bistable or multistable biochemical systems. Recently, parameter optimization procedures have been applied to calibrate oscillatory processes and to search for oscillatory and bistable massaction networks (Pitt and Banga, 2019; Porubsky and Sauro, 2019). However, for biological systems with multiple attractors – such as steadystate focus nodes, limit cycles, or more complex chaotic structures, which can coexist at the same parameter values – the use of standard parameter optimization pipelines that fit timeseries data might be extremely challenging. Also, systems of higher than 2D dimensions can potentially exhibit chaotic dynamics, typically happening according to one of the following scenarios: (i) perioddoubling bifurcations, (ii) a transition to chaos through intermittency, and (iii) chaotization of quasiperiodic dynamics (i.e., through the destruction of multidimensional tori) (Argyris et al., 1993; Kuznetsov, 2012).
Key aspects of the behavior of 2D dynamic systems are determined by the shapes of the nullclines (Tsyganov et al., 2012). Likewise, when a highdimensional system does not exhibit the chaotic behavior, the assembly of nullclines that are projections of quasisteadystate (QSS) trajectories into different 2D spaces help us understand the intricate system dynamics. To obtain the QSS trajectories, we change only one variable, whereas all other variables became implicit functions of that variable. Parameter optimization software packages developed for systems biology do not operate with nullclines. A critical feature of our computational pipeline is a combination of intensive sampling of the parameters space, identification of the number and stability of steady states, and the calculation and analysis of nullclines. Only when the key properties of nullclines are determined, the standard fitting software [e.g., BioNetFit (Thomas et al., 2016; Mitra et al., 2019) or PEPSSBI (Degasperi et al., 2017)] is used for finetuning of the parameters. After the dynamic regimes of the reduced 2D systems are established, the behavior of the original, highdimensional system is investigated in the vicinity of attractors found for reduced systems. This pipeline allows developing predictive models of such systems in a semiautomatic regime.
A schematic diagram highlighting the main steps of the proposed protocol is shown in Figure 1. The process consists of five successive stages:
1. Build an ODE model of the signaling network and constrain the parameter ranges using the available data on the protein abundances and kinetic constants.
2. Using a software package, e.g., the DYVIPAC package (Nguyen et al., 2015), which uses the libroadrunner library (Somogyi et al., 2015) for importing and solving ODE models, sample the parameter space, determine the number and stability of steady states for a given set of parameters, and analyze the shapes of nullclines by projecting the QSS trajectories into different 2D planes.
3. Scan 2D parametric planes and determine the types of bifurcations that occur in the system.
4. If for given dynamic regimes the system trajectories are measured, finetune the parameters to minimize deviations of simulated trajectories from the measured time series data using local search algorithms [e.g., simplex from BioNetFit package (Mitra et al., 2019)].
5. Based on the established dynamic regimes that often differ in distinct spatial regions, formulate a PDE system and perform spatial calculations.
Figure 1. Protocol overview. Consecutive steps of the protocol are indicated by arrows.
Equipment

Laptop computer
Proposed protocols can be implemented on a laptop computer; however, parameter sampling procedure is timeconsuming, and using a computational cluster is advantageous. For example, sampling parameter space of the model (BoladoCarrancio et al., 2020) to obtain a comprehensive list of different regimes took approximately 5 days at a workstation with 4core Intel^{®} Core^{TM} i7 Processor. The highquality scan of 2D parametric plane takes approximately 5 hours given the same computational power. Because the sampling procedures are very well parallelized, the computational time is inversely proportional to the number of available cores.
Since most computational clusters use Linuxbased operating systems, we further refer to software packages mainly developed for Linux.
Software

Python and python libraries, including matplotlib (Hunter, 2007), NumPy (Harris et al., 2020), and scipy (Virtanen et al., 2020)

Python version of DYVIPAC software package (Nguyen et al., 2015)

Software for parameter fitting of ODE systems, e.g., BioNetFit (Thomas et al., 2016; Mitra et al., 2019), PEPSSBI (Degasperi et al., 2017), etc.

Software for performing spatiotemporal calculations of reactiondiffusionconvection models, e.g., FiPy (Guyer et al., 2009), OpenFOAM (Jasak, 2009; OpenCFD, 2009) or Virtual Cell (Loew and Schaff, 2001)
Procedure

Build an ODE model of the signaling network based on the available topology of interactions between nodes. Implement the model in a software package that supports export to the SBML format, such as COPASI (Hoops et al., 2006) or BioNetGen (Blinov et al., 2004; Harris et al., 2016).

If the kinetics of interactions between two nodes is known, this network connection is modeled mechanistically.

For other connections, the interaction kinetics might not yet be exact mechanistically characterized, or these connections operate via intermediate species that are not included in the network under study. These connections are modeled using hyperbolic multipliers α_{ij} that specify the negative or positive influence of the activity (x_{j}) of node j on the activity (x_{i}) of node i (Tsyganov et al., 2012; BoladoCarrancio et al., 2020),
The coefficient γ_{ij}>1 indicates activation; γ_{ij}<1 inhibition; and γ_{ij}=1 denotes the absence of regulatory interactions, in which case the modifying multiplier α_{ij} equals 1. K_{ij} is the activation or inhibition constant.

Specify the parameter ranges using the data on the protein abundances and kinetic constants and assign numerical values to model parameters.

Nondimensionalize the model to reduce the number of parameters (Barenblatt, 2003).

Export the model to SBML format.

Determine dynamical regimes that can be observed in the system.

Import the sbmlfile of the model into the DYVIPAC package. Sample parameter space of the model using the DYVIPAC package, which determines the number and stability of steady states for a given set of parameters.

Rank protein abundances in a model in descending order and pick two proteins (nodes) with the highest abundance. Using the QSS approximation for the other protein activities (Sauro and Kholodenko, 2004; Eilertsen and Schnell, 2020), derive a 2D model from the initial ODE model.

Based on the specification of the steady states and the analysis of the shapes of nullclines (Figure 2, for example), determine the dynamic regimes and validate these regimes by the integration of ODEs.
Figure 2. Examples of the nullcline analysis. Nullclines and vector fields calculated for 2D ODE system are derived from a fivedimensional ODE system using a quasisteadystate approximation. Circles show stable steady states; triangles represent unstable steady states. Red and blue curves are nullclines for variables x_{1} and x_{2}, respectively. Green line represents trajectories of limit cycles projected from the original fivedimensional system to 2D space of x_{1} and x_{2} [see BoladoCarrancio et al. (2020) for details].

If more than two proteins have high and comparable abundances, make a list of pairs of these proteins and apply the above point c) to the ODE system for each pair.

Scan 2D parametric planes and determine the types of bifurcations that can occur in the system.

Scan the parameter planes and determine the borders that separate different dynamic regimes (see Figure 3, for example). This step can predict the previously unobserved dynamic regimes, which must be validated in the experiments.
Figure 3. Examples of 2D parametric scans. Different 2D parametric diagrams obtained using scanning of 2parameter planes. Different colors indicate different dynamical regimes. Black lines represent borders between these regimes where bifurcations happen (see BoladoCarrancio et al., 2020 for details and Figure 2 there).

The analysis of the changes in the number and stability of steady states can reveal local bifurcations, such as the AndronovHopf or saddlenode bifurcations (Kuznetsov, 2004). However, nonlocal bifurcations might exist in a 2D system where a limit cycle can appear or vanish (Nekorkin, 2015). Scan the parameters for dynamic regimes that contain focuses to determine the appearance or annihilation of limit cycles.

If there are timeseries data, finetune model parameters to reproduce the system trajectories for a biologically relevant regime. This step uses standard software packages for parameter optimization (Glover et al., 2000; Maiwald and Timmer, 2008; Kreutz and Timmer, 2009; Penas et al., 2015; Thomas et al., 2016; Degasperi et al., 2017; Mitra et al., 2018 and 2019).

Derive a PDE system by adding the terms describing mass transfer processes such as diffusion and convection to the ODE system studied (Tsyganov et al., 2012; BoladoCarrancio et al., 2020). Perform spatiotemporal simulations using specific software packages, such as OpenFOAM, FiPy, or Virtual Cell.
A set of examples that illustrate the analysis of the dynamics of the RhoARac1 signaling network model and python scripts can be found in “examples.zip” archive.
Acknowledgments
This protocol was adapted from Tsyganov et al. (2012) and Boladocarrancio et al. (2020). This work was supported by NIH/NCI grant R01CA244660 and EU grant NanoCommons (grant no. 731032).
Competing interests
The authors declare no competing interests.
References

Argyris, J., Faust, G., and Haase, M. (1993). Routes to Chaos and Turbulence. A Computational Introduction. In: Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences 344: 207234.

Barenblatt, G. I. (2003). Scaling. Cambridge University Press, Cambridge.

Blinov, M. L., Faeder, J. R., Goldstein, B. and Hlavacek, W. S. (2004). BioNetGen: software for rulebased modeling of signal transduction based on the interactions of molecular domains. Bioinformatics 20(17): 32893291.

BoladoCarrancio, A., Rukhlenko, O. S., Nikonova, E., Tsyganov, M. A., Wheeler, A., GarciaMunoz, A., Kolch, W., von Kriegsheim, A. and Kholodenko, B. N. (2020). Periodic propagating waves coordinate RhoGTPase network dynamics at the leading and trailing edges during cell migration. Elife 9: e58165.

Degasperi, A., Fey, D. and Kholodenko, B. N. (2017). Performance of objective functions and optimisation procedures for parameter estimation in system biology models. NPJ Syst Biol Appl 3: 20.

Eilertsen, J. and Schnell, S. (2020). The quasisteadystate approximations revisited: Timescales, small parameters, singularities, and normal forms in enzyme kinetics. Math Biosci 325: 108339.

Glover, F., Laguna, M. and Martí, R. (2000). Fundamentals of scatter search and path relinking. Control and Cybernetics 29: 652684.

Guyer, J.E., Wheeler, D. and Warren, J. A. (2009). FiPy: Partial Differential Equations with Python. Comput Sci Eng 11(3): 615.

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., Del Rio, J. F., Wiebe, M., Peterson, P., GerardMarchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C. and Oliphant, T. E. (2020). Array programming with NumPy. Nature 585(7825): 357362.

Harris, L. A., Hogg, J. S., Tapia, J. J., Sekar, J. A., Gupta, S., Korsunsky, I., Arora, A., Barua, D., Sheehan, R. P. and Faeder, J. R. (2016). BioNetGen 2.2: advances in rulebased modeling. Bioinformatics 32(21): 33663368.

Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., Xu, L., Mendes, P. and Kummer, U. (2006). COPASIa COmplex PAthway SImulator. Bioinformatics 22(24): 30673074.

Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Comput Sci Eng 9(3): 9095.

Jasak, H. (2009). OpenFOAM: Open source CFD in research and industry. Int J Nav Archit Ocean Eng 1(2): 8994.

Kreutz, C. and Timmer, J. (2009). Systems biology: experimental design. FEBS J 276(4): 923942.

Kuznetsov, Y. A. (2004). Elements of applied bifurcation theory. (3rd Edition). Springer.

Kuznetsov, S. P. (2012). Hyperbolic Chaos: A Physicist’s View. Springer.

Loew, L. M. and Schaff, J. C. (2001). The Virtual Cell: a software environment for computational cell biology. Trends Biotechnol 19(10): 401406.

Maiwald, T. and Timmer, J. (2008). Dynamical modeling and multiexperiment fitting with PottersWheel. Bioinformatics 24(18): 20372043.

Mitra, E.D., Dias, R., Posner, R.G., and Hlavacek, W.S. (2018). Using both qualitative and quantitative data in parameter identification for systems biology models. Nat Commun 9: 3901.

Mitra, E. D., Suderman, R., Colvin, J., Ionkov, A., Hu, A., Sauro, H. M., Posner, R. G. and Hlavacek, W. S. (2019). PyBioNetFit and the Biological Property Specification Language. iScience 19: 10121036.

Nekorkin, V. I. (2015). Chapter 9: Bifurcations of Limit Cycles. Saddle Homoclinic Bifurcation. In: Introduction to Nonlinear Oscillations. 113122. John Wiley & Sons, Inc.

Nguyen, L. K., Degasperi, A., Cotter, P. and Kholodenko, B. N. (2015). DYVIPAC: an integrated analysis and visualisation framework to probe multidimensional biological networks. Sci Rep 5: 12569.

OpenCFD. (2009). OpenFOAM. The Open Source CFD Toolbox. User Guide. OpenCFD Limited.

Penas, D. R., Banga, J. R., González, P. and Doallo, R. (2015). Enhanced parallel Differential Evolution algorithm for problems in computational systems biology. Appl Soft Comput 33: 8699.

Pitt, J. A. and Banga, J. R. (2019). Parameter estimation in models of biological oscillators: an automated regularised estimation approach. BMC Bioinformatics 20(1): 82.

Porubsky, V. L., and Sauro, H. M. (2019). Application of parameter optimization to search for oscillatory massaction networks using python. Processes 7(163): 163.

Sauro, H. M., and Kholodenko, B. N. (2004). Quantitative analysis of signaling networks. Prog Biophys Mol Biol 86(1): 543.

Somogyi, E. T., Bouteiller, J. M., Glazier, J. A., Konig, M., Medley, J. K., Swat, M. H. and Sauro, H. M. (2015). libRoadRunner: a high performance SBML simulation and analysis library. Bioinformatics 31(20): 33153321.

Thomas, B. R., Chylek, L. A., Colvin, J., Sirimulla, S., Clayton, A. H., Hlavacek, W. S. and Posner, R. G. (2016). BioNetFit: a fitting tool compatible with BioNetGen, NFsim and distributed computing environments. Bioinformatics 32(5): 798800.

Tsyganov, M. A., Kolch, W. and Kholodenko, B. N. (2012). The topology design principles that determine the spatiotemporal dynamics of Gprotein cascades. Mol Biosyst 8(3): 730743.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P. and SciPy, C. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17(3): 261272.
Copyright
Rukhlenko and Kholodenko. This article is distributed under the terms of the
Creative Commons Attribution License (CC BY 4.0).
How to cite: Readers should cite both the Bioprotocol article and the original research article where this protocol was used:
 Rukhlenko, O. S. and Kholodenko, B. N. (2021). Modeling the Nonlinear Dynamics of Intracellular Signaling Networks. Bioprotocol 11(14): e4089. DOI: 10.21769/BioProtoc.4089.
 BoladoCarrancio, A., Rukhlenko, O. S., Nikonova, E., Tsyganov, M. A., Wheeler, A., GarciaMunoz, A., Kolch, W., von Kriegsheim, A. and Kholodenko, B. N. (2020). Periodic propagating waves coordinate RhoGTPase network dynamics at the leading and trailing edges during cell migration. Elife 9: e58165.