Modeling the Nonlinear Dynamics of Intracellular Signaling Networks.

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 fine-tune 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 quasi-steady-state 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. ©Copyright Rukhlenko et al.

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 ; Bolado- Carrancio ). A variable 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 ( or , respectively), the element 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 ; Maiwald and Timmer, 2008; Kreutz and Timmer, 2009; Penas ; Thomas ; Degasperi ; Mitra 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 mass-action networks (Pitt and Banga, 2019; Porubsky and Sauro, 2019). However, for biological systems with multiple attractors – such as steady-state focus nodes, limit cycles, or more complex chaotic structures, which can co-exist at the same parameter values – the use of standard parameter optimization pipelines that fit time-series 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) period-doubling bifurcations, (ii) a transition to chaos through intermittency, and (iii) chaotization of quasi-periodic dynamics (i.e., through the destruction of multidimensional tori) ( Argyris ; Kuznetsov, 2012).

Key aspects of the behavior of 2D dynamic systems are determined by the shapes of the nullclines ( Tsyganov ). Likewise, when a high-dimensional system does not exhibit the chaotic behavior, the assembly of nullclines that are projections of quasi-steady-state (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 ; Mitra ) or PEPSSBI ( Degasperi )] is used for fine-tuning of the parameters. After the dynamic regimes of the reduced 2D systems are established, the behavior of the original, high-dimensional system is investigated in the vicinity of attractors found for reduced systems. This pipeline allows developing predictive models of such systems in a semi-automatic regime.

A schematic diagram highlighting the main steps of the proposed protocol is shown in Figure 1. The process consists of five successive stages:

Figure 1.

Protocol overview.

Consecutive steps of the protocol are indicated by arrows.

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 ), which uses the libroadrunner library ( Somogyi ) 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, fine-tune 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 )].

5. Based on the established dynamic regimes that often differ in distinct spatial regions, formulate a PDE system and perform spatial calculations.

Laptop computer

Proposed protocols can be implemented on a laptop computer; however, parameter sampling procedure is time-consuming, and using a computational cluster is advantageous. For example, sampling parameter space of the model (Bolado- Carrancio ) to obtain a comprehensive list of different regimes took approximately 5 days at a workstation with 4-core Intel® CoreTM i7 Processor. The high-quality 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 Linux-based operating systems, we further refer to software packages mainly developed for Linux.

Python and python libraries, including matplotlib (Hunter, 2007), NumPy ( Harris ), and scipy ( Virtanen )

Python version of DYVIPAC software package ( Nguyen )

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

Software for performing spatiotemporal calculations of reaction-diffusion-convection models, e.g., FiPy ( Guyer ), OpenFOAM (Jasak, 2009; OpenCFD, 2009) or Virtual Cell (Loew and Schaff, 2001)

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 ) or BioNetGen ( Blinov ; Harris ).

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 that specify the negative or positive influence of the activity () of node j on the activity () of node i ( Tsyganov ; Bolado- Carrancio ),

The coefficient indicates activation; inhibition; and denotes the absence of regulatory interactions, in which case the modifying multiplier equals 1. 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.

Non-dimensionalize 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 sbml-file 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 five-dimensional ODE system using a quasi-steady-state approximation. Circles show stable steady states; triangles represent unstable steady states. Red and blue curves are nullclines for variables x1 and x2, respectively. Green line represents trajectories of limit cycles projected from the original five-dimensional system to 2D space of x1 and x2 [see Bolado-Carrancio 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 2-parameter planes. Different colors indicate different dynamical regimes. Black lines represent borders between these regimes where bifurcations happen (see Bolado- Carrancio 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 Andronov-Hopf or saddle-node bifurcations (Kuznetsov, 2004). However, non-local 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 time-series data, fine-tune model parameters to reproduce the system trajectories for a biologically relevant regime. This step uses standard software packages for parameter optimization ( Glover ; Maiwald and Timmer, 2008; Kreutz and Timmer, 2009; Penas ; Thomas ; Degasperi ; Mitra 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 ; Bolado- Carrancio ). Perform spatiotemporal simulations using specific software packages, such as OpenFOAM, FiPy, or Virtual Cell.

Procedure

A set of examples that illustrate the analysis of the dynamics of the RhoA-Rac1 signaling network model and python scripts can be found in “examples.zip” archive.

dx_idt=f_ix_1,…,x_n,   i=1, 2, …, n

(1).

α_ij=1+γ_ijx_j/K_j1+x_j/K_j

(2).