Steady-state kinetic modeling constrains cellular resting states and dynamic behavior

Abstract

A defining characteristic of living cells is the ability to respond dynamically to external stimuli while maintaining homeostasis under resting conditions. Capturing both of these features in a single kinetic model is difficult because the model must be able to reproduce both behaviors using the same set of molecular components. Here, we show how combining small, well-defined steady-state networks provides an efficient means of constructing large-scale kinetic models that exhibit realistic resting and dynamic behaviors. By requiring each kinetic module to be homeostatic (at steady state under resting conditions), the method proceeds by (i) computing steady-state solutions to a system of ordinary differential equations for each module, (ii) applying principal component analysis to each set of solutions to capture the steady-state solution space of each module network, and (iii) combining optimal search directions from all modules to form a global steady-state space that is searched for accurate simulation of the time-dependent behavior of the whole system upon perturbation. Importantly, this stepwise approach retains the nonlinear rate expressions that govern each reaction in the system and enforces constraints on the range of allowable concentration states for the full-scale model. These constraints not only reduce the computational cost of fitting experimental time-series data but can also provide insight into limitations on system concentrations and architecture. To demonstrate application of the method, we show how small kinetic perturbations in a modular model of platelet P2Y(1) signaling can cause widespread compensatory effects on cellular resting states.

Figure 1. Example structure of model topology, ODEs, and concentration space.

(A) The model topology defines the state transitions (arrows) and rate equations (f) that determine how molecules are interconverted. This example model is organized into three overlapping modules, with molecules and each occurring in two modules. Corresponding (B) ODEs and (C) concentration space for the example topology in panel A. Each of the 7 molecules occupies a separate linear dimension, with each module comprising a subspace of the full 7-dimensional space. Modules that share a common molecule have intersecting subspaces.

Figure 2. Steps in dimensionality reduction of steady-state modules and example from platelet signaling model.

(A) Steps in dimensionality reduction of kinetic modules: (1) Restrict value ranges for each c to physiologically realistic ranges. (2) Compute multiple steady-state solutions to the model ODEs using initial guesses sampled randomly from the defined distribution. (3) Reduce the dimensionality of the steady-state solution set by PCA. (B) Results obtained from modular reduction 4 kinetic modules in a platelet signaling model. For each module, a fixed topology was combined with initial guesses from the defined distribution and simulated until equilibrium was reached () using 10⁹ initial guesses for . Specific concentrations within these steady-state solutions were compared to experimentally measured values, and solutions with low error (±10% of known concentration values) for these elements were selected as “points” in the steady-state concentration space. PCA was then applied to transform these points to a new coordinate set that maximally covers the space of steady-state solutions.

Figure 3. Assembly of full model from steady-state modules.

(A) The full model is assembled by combining PCA-reduced, steady-state solution spaces from each module into a combined steady-state solution space. This global space is searched for full-length, steady-state solution vectors that satisfy both the steady-state requirements of each module and the desired time-dependent properties when the steady-state is perturbed (in this example, by increasing the concentration of the signaling molecule ADP and measuring the change in intracellular Ca²⁺ concentration). A simple linear constraint is imposed for every pair of modules that share a common molecule to ensure that steady-state solutions are consistent. (B) To assemble the platelet signaling model, a set of 16 PC vectors representing all 72 unknown variables in the model were used as search directions in a global optimization routine. The global solution space was searched for models with accurate dynamic behavior using experimental time-series data for ADP-stimulated Ca²⁺ release. Species are grouped according to compartment. Color values correspond to molar concentrations (mol/L or mol/m²) or as indicated: *DTS species (mol L⁻¹). †Extracellular species (mol L⁻¹). ‡DTS volume (L). §PM leak conductance/area (S m⁻²).

Figure 4. Shifts in steady-state profiles caused by kinetic perturbations.

The steady-state platelet model was perturbed by changing selected kinetic parameters (±10-fold) and simulating for 1 h (left panels). After approaching a new steady state, the model concentrations and fluxes were determined relative to their original steady-state values and colored according to fold-change (right panels). Green indicates no change (NC) relative to initial flux/concentration. Red indicates a relative increase and blue indicates a relative decrease. Note that the color scale in each panel is normalized separately to maximize distinctions in fold change. New steady states were achieved after a (A) 10-fold increase in Ca²⁺ release through open IP₃R channels ([28]), (B) 10-fold decrease in PKC-mediated inhibition of PLC-β, and (C) 10-fold increase in PIP₂ hydrolysis (10-fold increase in k _cat of hydrolysis). Reactions with perturbed rate constants are circled and correspond to reaction mechanisms from [ref. 14]. (A) Ca²⁺ _dts → Ca²⁺ _i, (B) PKC*+PLC-β → PKC*+pPKC-β, (C) PLC-β*+PIP₂ → PLC-β*+IP₃+DAG.