Preprint

Article

Linear and Nonlinear Modes & Data Signatures in Dynamic Systems Biology Models

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

Joseph DiStefano^*

This version is not peer-reviewed

Submitted:

26 July 2023

Posted:

28 July 2023

You are already at the latest version

Alerts

Abstract

The particulars of stimulus-response experiments done on dynamic biosystems clearly limit what one can learn and validate about their structural interconnectivity (topology), even when collected kinetic output data are perfect (noise-free). As always, available access ports and other data limitations rule. For linear systems, exponential modes, visible and hidden, play an important role in understanding data limitations, embodied in what we call dynamical signatures in the data. We show here how to circumscribe and analyze modal response data in compartmentalizing model structures – so that modal analysis can be used constructively in systems biology model building – for nonlinear (NL) as well as linear biosystems. We do this by developing and exploiting the modal basis for dynamical signatures in hypothetical (perfect) input-output data associated with a structural model – one that includes inputs and outputs explicitly – and for NL as well as linear biosystems. The methodology establishes model dimensionality (size, complexity) from particular data sets; helps select among multiple candidate models (model distinguishability); helps in designing new experiments to extract “hidden” structure; and helps to simplify (reduce) models to their essentials. For NL biosystems, results are not as comprehensive, similarly informative about their dominant dynamical properties, and unified with linear models on invariant 2-dimensional manifolds in phase space. Some automation of these highly technical aspects of biomodeling is also introduced.

Keywords:

Subject: Engineering - Bioengineering

Introduction

Over the last few decades, interest and progress in mechanistic math modeling has burgeoned in the realm of dynamic systems biology. The subject literature has grown substantially, with much greater attention being paid to quantitative modeling methodology and applications, particularly in the university life science and pharmaceutical communities. Hindsight suggests that this focused diligence was to be expected. Development of vastly superior molecular and cell biology measurement tools and the numerous biological discoveries emanating from the explosive growth of data provided by these tools was undoubtedly a major if not the impetus. The common understanding is that explicit quantitative models transform data into useful constructs, ideas and discoveries. And the literature is growing with papers on modeling methodologies specialized for molecular and cellular level modeling and textbooks that are consolidating them e.g., (Palsson 2006, Alon 2007, Klipp 2009, Voit 2012, DiStefano III 2014, Distefano 2019).

Multicompartmental modeling – despite being prominent in the quantitative physiology and pharmacology literature for over half a century – remains an area underappreciated by systems biology modelers. Indeed, the subject of compartmental modeling occupies few if any pages in modern texts on systems biology, as for most of the references above. An overarching goal in this paper is, first, to make the case that this oversight may be hampering progress in the field – at least in the sense that readily available and effective tools for mechanistic modeling are being neglected. To accomplish this, I show the ease with which dynamic systems input-output (I-O) data can be exploited in structural, mechanistic model building, first using linear systems analysis modal concepts merged with linear compartmental modeling. Then I extend and unify these notions to nonlinear dynamic systems biology modeling, by merging and extending concepts in structural mechanics with and to biochemical kinetics. First, some theory and examples.

Linear Modes

Modal concepts are classical in linear system theory, e.g see (Chen 1985). Consider the general class of linear time-invariant ordinary differential equation (ODE) models, with inputs and outputs specified integral to the model M, as in Figure 1 and Eq. (1):

\dot{q} (t) = K q (t) + B u (t), y (t) = C q (t), t \geq 0

(1)

Symbols q, u and y are n-, r- and m-vectors of time-varying state variables (system), inputs and outputs. K, B and C are n by n system, n by r input and m by n output matrices of constants. Multicompartmental (MC) models are a special case, with

K \equiv [k_{i j}]

having the properties: k_ij ≥ 0 for i ≠ j and

\sum_{i = 0}^{n} k_{i j} = 0, j = 1, ..., n

(Rubinow 1975, Godfrey 1983).

The component terms of the solution

y (t) = C q (t) = C e^{K t} q (0)

of the homogeneous equation

\dot{q} = K q

, i.e., the response to initial conditions (ICs)

q (0)

(input u ≡ 0) are the modes of model (1). This fundamental set of solution terms (modes) are n linearly independent exponential functions of the form

t^{k} e^{λ_{j} t}

(k = 0, 1,…, n-1 and j = 1, 2,…, n), each satisfying

\dot{q} = K q

. The λ_j are eigenvalues of K. If they are distinct, the n modes are:

y_{1} = e^{λ_{1} t}, \dots, y_{n} = e^{λ_{n} t}

. If they are repeated, then for each λ_i of multiplicity n_i (i.e., n_i equal to λ_i), there are n_i modes:

e^{λ_{i} t}, t e^{λ_{i} t}, \dots, t^{n_{i} - 1} e^{λ_{i} t},

also totaling n.

Importantly, n is the minimum dynamical dimension of model (1), i.e., the model cannot be reduced to a smaller number of ODEs that completely characterize dynamical responses

y (t) = C q (t) = C e^{K t} q (0)

. This is paramount in tying model with data, as we see below. The eigenvalues λ_j can be real or complex numbers. When they are complex, the modes are oscillatory.

We now give all this a geometric interpretation, to tie it to discussion of nonlinear (NL) models and modes later, Each mode actually has a shape, represented by an eigenvector component v_j, associated with the frequency (eigenvalue) component λ_j, so it can be represented geometrically as a two-dimensional function

c_{j} v_{j} t^{k} e^{λ_{j} t}

each with its own eigenvalue λ_j and (n-dimensional) eigenvector v_j. The intrinsic parameters of

y (t) = C q (t) = C e^{K t} q (0)

are then the eigenvalues and eigenvectors of K, i.e., a set of 2n structural invariants of K. Structural invariants are one of many smallest sets of parameters completely determining or characterizing model dynamics (Chen 1970).

Example 1: Simplest Case – Distinct Eigenvalues

The dynamical shape of the solution modes, for recognition purposes as well as math analysis, depends in part on relationships among the n eigenvalues λ₁, λ₂, …, λ_n of system matrix K. If they are distinct (none are repeated) and the corresponding eigenvectors are v₁, v₂, …,v_n, the general solution of

\dot{q} = K q

can be written:

q (t) = \sum_{i = 1}^{n} c_{i} v_{i} e^{λ_{i} t}

(2)

The c_i are chosen to match the initial conditions (ICs) on q(t). For example, if only

q_{1} (0)

≠ 0 and all others ICs are zero, then Eq. (2) gives:

q (t) = {[q_{1} (0) \sum_{i = 1}^{n} v_{1} e^{λ_{i} t} 0 \dots 0]}^{T}

. If the eigenvalues are not distinct, the solution is more complex, but still multiexponential, and each repeated eigenvalue has its own mode, e.g.,

c_{1} e^{λ t}, c_{2} t e^{λ t}

are two different modes for repeated eigenvalue λ. MC models can have repeated and/or complex eigenvalues in their responses.

We can generalize this for nonhomogeneous system inputs u(t) ≠ 0 in Eq. (1). Since the exponential mode response to ICs is the same as response to impulse inputs, i.e.,

B u (t) = q (0) δ (t)

in Eq. (1), the system response to step, ramp, exponential, etc (different) inputs – by linearity – also will be a sum of exponential modes and additional terms inherited from u(t).

A most important feature of linear modes is that they are a minimal set of invariants of the dynamical system and, for any set of ICs or input(s), the modes accurately capture essential system dynamics. Among other things, this means linear ODEs of order higher than n can be accurately reduced (simplified) to order n if the ODE has only n modes, and this simplified model can be represented by a set of n modes, each represented geometrically as a 2-dimensional function of the eigenvalue and eigenvector associated with each model. Whether all these modes are visible or hidden in data is the question we address next – the key to unraveling structure from data and designing new input-output (I-O) experiments that reveal more.

Modes as Data Signatures

As implied above, both visible and hidden modes play important roles in understanding dynamical signatures in data, and thus the structure of systems from which they emanate. Using a possibly familiar analogy with large-scale data analysis in bioinformatics or computational statistics, finding the modes in dynamic system data is akin to finding: the principal components (principal component analysis – PCA), or singular values (singular-value decomposition – SVD), of a large data matrix. This means visualizing a multivariate dataset – typically expressed in a high-dimensional data space, from its “most informative” viewpoint – its principal components – in a lower-dimensional space. Distinct modes of chords and scales in music are another albeit less mathematical example, with each mode sounding different but having a different arrangement of the same notes played in particular orders.

This is important because we can use these notions for developing dynamic system model structure and establishing the dimensionality of the model observable in data, in terms of modes associated with n exponential solutions in the data. In brief, visible modes in data collected from a system are visibly associated or identified with. visible state variables in its model. Hidden modes in the data are associated with hidden state variables[1], i.e., those that have no signature in the data. They may be in the system, but they are not visible in the data. Different experiments (input-output configurations) are needed to discover them in new data.

To tie these relationships to model structure, in particular multicompartmental structural models, we need some definitions.

State Variables, Compartments, Directed, Graphs, Pools & Species

We formalize and further clarify our primary nomenclature, paraphrasing classical definitions collectively, in the contextual domain of biological systems, e.g., see (Rubinow 1975, Godfrey 1983, Jacquez 1996). [I found this section more difficult than other sections]

A compartment is as an amount of material X which acts kinetically (within a dynamic system) in a homogeneously distinct manner. Homogeneously distinct are the key words here. This means entities within a compartment are indistinguishable from each other, they are “mixed completely”. If they are not, they comprise more than one compartment. A compartment is open if it “leaks” to the environment; it is closed if it does not (Figure 2). Most (not all) compartmental models of real systems have at least one open compartment.

Every (physical) space containing compartments has a presumed volume in which the material substance (e.g., chemical “species”, etc) is homogeneously distributed. This compartment volume V_i might be measurably real, or it might be virtual – called an equivalent volume.

In Other Words, in chemical, biological or other physical spaces, a compartment is an idealized store of a substance, characterized by its environment, physiochemical properties, or both. All the substance in a given form qualifies as a compartment, e.g., active or inactive, or receptor-bound or unbound form – one compartment for each – because they generally have different dynamics. Alternatively, all of the substance in a homogeneously distributed (well-mixed) location, or all of it in a given form and location, each qualify as compartments, e.g., the mass of unbound drug D_F (or hormone, or lipid, or ion) in blood, or the mass of receptor-bound drug D_R in blood, or either in liver (or other organ, or organelle). Similarly, the total drug mass

D_{F} + D_{R}

in blood or liver (or other organ, or organelle), are each valid compartments in principle. Similarly, in modeling epidemiological systems, Susceptible, Exposed, Infected and Recovered (SEIR) populations are typically each considered valid, homogeneous compartments.

A multicompartmental (MC) model consists of a directed graph of two or more compartments interconnected so there can be exchange of material, or change of state (into two or more compartments), among some or all of them, and possible leaks from any of them (i.e., Figure 2). Exchange may occur by material transgressing some physical barrier (e.g., by diffusion or membrane transport), or by undergoing some physical or chemical transformation (e.g., metabolism, inactive to active state, etc). Each compartment is represented by a single state variable in a multicompartmental (MC) model. Changes in location of material generate different compartments (different state variables) for the same material, typically with different volumes. Changes in material form (e.g., precursor to product, carbohydrate to metabolite, etc) generate different compartments (different state variables) – possibly in the same “space” or “container” and with the same volume, and so on.

Circles are commonly used to represent compartments, arrows unidirectional transfers – which include inputs, influxes and effluxes (leaks) – and dashed lines, with a small open-circle endpoint, denote measured outputs (Figure 3). Directed graphs (digraphs) also are used to represent MC models, with nodes (vertices) representing compartments and edges as the arrows interconnecting compartments (nodes) or the environment, as in Fig. 3. For linear systems, arrows and edges are typically represented by constant parameters (e.g., fractional transfer rates k_ij). For nonlinear systems, arrows and edges designate functions of constant parameters (e.g., A, B) and state variables q_i (e.g., saturable functions, like A/(B + q_i) (Distefano 2019).

A pool is the total amount of a substance in a system or subsystem (Bergner 1967), defined purposely loosely to suit the variety of ways different investigators were using it. We are a bit more precise here, defining pool – similarly, as a mass or volume size measure for a real or abstract space or material, emphasizing that it is not necessarily uniformly (homogeneously) distributed (Rescigno 1990). A compartment may be a pool, but not necessarily vice versa. On the other hand, we can speak of the “pool of X in compartments i and j” or the “pool of X in the whole organism” (or “whole compartmental model or system”) because the term pool circumvents the homogeneity requirement on any of these entities. Pool Models are particularly useful for steady state flux balance analysis (FBA) (DiStefano 2013).

Species is a term commonly used in chemistry and molecular biology, with nomenclature embodied in the latest version (Level 3) of the Systems Biology Markup Language SBML, a “bridge language (lingua franca) for communication or exchange of mathematical models across programming tools” (sbml.org). Chemical species are defined as distinct molecular entities sharing the same chemical properties, and same chemical properties effectively means homogeneous. In a chemical or biological system kinetic (dynamical) model, distinct species are thus represented as distinct state variables – the classical definition of (homogeneous) compartment in a dynamic system model[2].

A multicompartmental (MC) model can have any number of compartments. But its observable (visible) dimensionality is limited by the locations of inputs and outputs and thus the particular (noise-free) data available. This means the granularity – or extent – of multicompartmentalization is experiment dependent. The minimum dimension – the smallest number of possible compartments – is equal to the number of its distinct modes visible in its output reponse – for any particular specification of input and output locations.

Modes, Compartments and Data: Modes in Data = Minimum Compartment Number

How do we characterize (or build) system dynamical structure from temporal kinetic input-output (I-O) data collected from accessible I-O ports? We have seen above that the modes of this data are the fundamental solutions of the ODEs representing the model, one for each of the state variables; and any zero-input response – i.e., response due only to ICs – is a linear combination of these modes. Thus, for linear multicompartmental models, the number of modes is equal to the number of state variables, which means mode-count ≡ minimum compartment-count.

The actual number of system compartments – identified as such or not – is typically greater than the minimum needed to represent the dynamics, in many cases very much larger than this minimum number. Probed compartment(s) are likely included in the count, but – at least in principle – additional “hidden” and unprobed compartments need to be discovered – typically by probing different I-O ports – to establish actual dynamical dimensionality, as illustrated below. This can be challenging for the modeler, no matter what the goals are.

Finding Modes (Compartments) Visible in Output Data by Graphical Inspection

Not all modes and compartments (state variables) they represent are visible in a particular output y(t), in response to a particular input u(t) – i.e., an input-output port pair. But they may be visible from other I-O data ports. We take advantage of the fact that modes are intrinsically responses to initial conditions (ICs) or, equivalently, impulse inputs. The procedure is to systematically test each I-O port pair of a candidate model structure, counting the number of modes in each designated output compartment in response to an impulse input into the designated input compartment. Checking all possible compartment I-O pairs in this manner then gives all possible modes in the system, i.e., all possible compartments in the candidate model. This, incidentally, reveals new experiment designs for finding otherwise hidden compartments.

Analysis can be done mathematically, but it’s easier by systematic inspection of the compartmental model graph. This graphical approach exposes visually why full dimensionality of a model and the system it represents is difficult to establish. In the process, it also provides a more facile tool for discovering hidden structure in underlying system. The following two simple examples illustrate the procedure step-by-step. Additional, more complex examples can be found in (DiStefano 2013).

Example 2: 2-Compartment Candidate Model with 1-Input & 2-Outputs

In the simple model of Figure 4, we assume zero ICs and a single impulse input. We relax these assumptions, later, to discover how ICs and other inputs affect results. The k_ijs are nonnegative rate constants (time^-1) and V_is are volumes. We are interested in the number of modes visible in outputs 1 and 2 under various conditions. We evaluate this number first mathematically, then more simply using graphical-inspection.

The Math Solution: It can be shown (DiStefano 2013) that both outputs are a sum of two real and distinct exponential modes:

y_{1} (t) = \frac{1}{V_{1} (λ_{1} - λ_{2})} [(k_{02} + k_{12} + λ_{1}) e^{λ_{1} t} - (k_{02} + k_{12} + λ_{2}) e^{λ_{2} t}], y_{2} (t) = \frac{k_{21}}{V_{1} (λ_{1} - λ_{2})} (e^{λ_{1} t} - e^{λ_{2} t})

and the two exponents (eigenvalues) are:

λ_{1, 2} = (- (k_{02} + k_{12} + k_{21}) \pm \sqrt{{(k_{02} + k_{12} + k_{21})}^{2} - 4 k_{02} k_{21}}) / 2 \leq 0 for all k_{i j} \geq 0

We have exactly two distinct (exponential) modes measurable in each compartment for an impulse input for all k_ij > 0. Now suppose k₂₁ = 0 (and ICs remain zero), then λ₂ = − (k₀₂+k₁₂), which nulls λ₁ and the second exponential term in the equation for y₁, i.e., the solution for y₁ reduces to only one mode, and that for y₂ reduces to zero modes because the coefficient goes to zero (by L’Hopital’s rule). Thus, the model with k₂₁ = 0 will exhibit one mode in 1 and none in 2. That’s what the math tells us.

The Simpler Graphical Solution: We redo the problem by reasoning physically about each case, k₂₁ nonzero and k₂₁ zero, i.e., with and without the forward connection from compartment 1 to 2. The principle here is that each compartment, when it is stimulated by a signal entering it, elicits a distinct mode associated with that compartment in its dynamic response. Whether this effect is visible in an output depends on the locations of the I-O probes, and the compartment interconnections.

An input signal entering compartment 1 (an impulse input in this example) will first stimulate the first mode represented by the state variable – dynamic system compartment 1 itself. Then, if k₂₁ is nonzero, the signal travels to 2 via k₂₁, where it stimulates the second distinct mode, represented by the state variable – dynamic system compartment 2. Finally, with k₁₂ nonzero, material travels back to 1 via k₁₂, with a dynamic reflecting both the stimulus from 1 along k₂₁ and the corresponding dynamic due to stimulation of the distinct mode in 2. So we see two modes in both compartments. They go back and forth, feeding their signals forward and back, because they are connected reversibly.

Other Parameterizations: Now, using the same reasoning and inspection of the graph, we see that if k₁₂ = 0 and k₂₁ is nonzero, we find only one mode in compartment 1 and 2 in compartment 2, for the same input and ICs, because the mode associated with compartment 2 is not visible in 1, because of the lack of connection back to 1 (k₁₂ = 0). We discovered this using mathematics, above, but this graphical solution is easier.

Nonzero ICs: Going one step further, if we started with zero input, k₁₂ ≠ 0, k₂₁ ≠ 0, and nonzero ICs in both compartments, both would have two modes visible in them, because the ICs are equivalent to impulse inputs. And, unlike the zero IC case above, when k₂₁ = 0, compartment 1 would still show 2 modes, because mode 2 would be stimulated by the nonzero IC in 2 and travel to 1 along k₁₂.

So nonzero ICs can be as important as inputs in establishing model complexity, because – in essence – they are equivalent to impulse inputs and generate the same modal responses.

Example 3: 3- Compartment Candidate Model

We add a third compartment, the only one measured, and evaluate the number of modes visible in the output of compartment 3 in Figure 5, using only graphical reasoning. We again assume zero ICs everywhere, an impulse input and all k_ij > 0. In the previous example we established that two modes would be visible in compartment 2, if it were measured. Despite the fact that compartment 2 is not measured in this problem, there are still two “excited” modes in 2, and they travel to compartment 3 along nonzero path k₃₂. So, at least two modes are visible in compartment 3 measurements. But, compartment 3 is itself a dynamic state variable compartment, and when it is stimulated it generates its own mode. Stimulated by two modes coming from 2, it therefore exhibits three modes in its output. No math needed!

Different Inputs: Going one step further, if the input u₁ were an exponential function distinct from the fundamental modes of the model, e.g., a single exponential input, then four modes would be visible at y₃, because the input “mode” would also be visible in the output. Inputs are external to the model; therefore exponential inputs can be thought of as the output of a compartment external to the model, another mode! And so on….

Remarks: In principle, the results above do not depend on parameter values, only structure, assuming all k_ij > 0 in the models. In reality, the relative magnitudes of the parameter values matter, because – depending on parameter values, some modes may be very small and difficult to detect in real data. This is a problem for model parameter sensitivity analysis, beyond our scope here.

Impulse inputs (approximate brief-pulse inputs in practice) are in principle maximally informative for modal analysis. If the input were instead a step, or ramp, or any integral of the impulse, the modal responses would be qualitatively the same – still exponential, with the same frequency λ_i (eigenvalues), but with different shapes (eigenvectors). Arbitrary inputs also elicit the same albeit less obvious modal responses; and they would be more difficult to unravel. Nevertheless, they too can be accommodated if needed, by suitable (e.g., Fourier transform) analysis – also beyond our scope in this paper.

Automated Mode Detection Using Graph Theory Algorithms

The graphical mode-finding algorithm of the last section is readily automated using computer science graph theory tools. Either Breadth-First Search (BFS) or Depth-First Search (DFS) (Cormen 2009) can be used to count compartments with a directed path between any measurable compartments, i.e., those reachable from perturbed compartments, i.e., those receiving inputs or having nonzero initial conditions (ICs are the same as equivalent impulse inputs). Directed connectedness and the permutations of paths between inputs to different compartments and follow-up compartments is key to counting modes. Pseudo code[3] for the imbedded mode-counting algorithm of the last section is given in the APPENDIX. Notation G(V,E) means a graph with V vertices (nodes, compartments) and E edges (rate constant parameters).

Hidden Modes & Model Simplification

Hidden modes and compartments[4] also have an upside, in the context of model simplification. For state variable models, model reduction (aggregation) means finding a simpler representation (with fewer state variables) appropriate or useful for some purpose. Transfer functions (TFs) H_ij(s) between pairs of inputs and outputs (u_j and y_i) of multicompartmental models

\dot{q} = K q + B u, y = C q

are particularly convenient for this, because pole-zero cancellations in H_ij(s) reduce model dimensionality, and therefore state variable and compartment number between these I-O pairs. This means that a multicompartmental model has hidden compartments if its TF has at least one pole-zero cancellation. Conversely, it has no hidden compartments if there are no pole-zero cancellations (Fagarasan and DiStefano III 1986).

Example 4: 3- Compartment Candidate Model Transfer Function & Simplifications

The TF for the model of Example 3, Fig. 5, is computed by Laplace transform analysis as:

H_{31} (s) = \frac{k_{21} k_{32}}{s (s^{2} + (k_{21} + k_{12} + k_{02} + k_{32}) s + k_{21} (k_{02} + k_{32}))}

(3)

No pole-zero cancellation is possible in (3) with all k_ij > 0, so model dimensionality cannot be reduced from 3. This is the case even if the turnover rates of compartments 1 and 2 are equal, i.e.,

{H_{31} (s) |}_{k_{21} = k_{12} + k_{02} + k_{32}} = \frac{k_{21} k_{32}}{s (s + k_{21} \pm \sqrt{k_{12} k_{21}})}

(4)

If both the input and output are instead in compartment 2 here (TF H₂₂), the TF is:

H_{22} (s) = \frac{s (s + (k_{12} + k_{02} + k_{32}))}{s (s^{2} + (k_{21} + k_{12} + k_{02} + k_{32}) s + k_{21} (k_{02} + k_{32}))}

(5)

The zero and pole at the origin in (5) cancel each other and the model is reduced to 2nd order. And if turnover rates in compartments 1 and 2 are equal, we still have a 2-dimensional model:

{H_{22} (s) |}_{k_{21} = k_{12} + k_{02} + k_{32}} = \frac{s + k_{21}}{{(s + k_{21})}^{2} - k_{21} k_{12}}

(6)

Although not the intention of this example, the number (two) of visible compartments in this H₂₂ example is available simply by inspection of the graph: the mode generated in compartment 3 cannot be returned to 2, because there is no path from 3 to 2.

Finally, H₂₂ can be reduced further to 1-dimension, by another pole-zero cancellation, for equal turnover rates of compartments 1 and 3 in Fig. 4, because their modes then have the same eigenvalues (frequencies) and thus the compartments cannot be distinguished from measurements in compartment 2. The math is:

{H_{11} (s) |}_{k_{21} = k_{12} + k_{02} + k_{32}} = \frac{s + k_{21}}{s + k_{21} \pm \sqrt{k_{12} k_{21}}}

(7)

Nonlinear Modes in Nonlinear Models

Nonlinear (NL) ODE models also have “transfer functions” (TFs), i.e., input-output solution response functions. Unlike linear models, however, their dynamic responses are dependent on inputs – input shapes or magnitudes. So we can’t write them only in terms of model parameters, like we do for linear model TFs. Nevertheless, the principle result above carries over to NL models. “Modes” in NL systems can cancel each other – at least approximately – or be vanishingly small in magnitude, or present for vanishingly small transient times – thereby rendering complex structure less visible in data characterized by specific inputs and outputs. Similarly simplified NL models, also with input-dependent dynamics, are characteristic of enzyme kinetic models. We tie the linear mode concept together here with a novel view of quasi-steady state approximations (QSSA) in nonlinear biochemical kinetics models – viewed as nonlinear modes with the same kinds of geometric properties as linear modes. The goal is to better understand NL systems and facilitate NL model building and analysis.

The field of mechanics has a developed theory of nonlinear modes associated with certain oscillatory phenomena in mechanical systems. Biological systems share many of the same dynamical properties, not recognized as such in the field. We make the connection precise here.

As noted earlier, linear modes are a minimal set of n-dimensional invariants of the dynamical system and, for any set of ICs or input(s), the modes accurately capture essential system dynamics. Among other things, this means linear ODEs of order higher than n can be accurately reduced (simplified) to order n, and this simplified model can be represented by a set of n two-dimensional modes in a state vector space. Mode shape is represented by eigenvector component v_j, associated with frequency (eigenvalue) component λ_j. The intrinsic parameters of measurement model

y (t) = C q (t) = C e^{K t} q (0)

for the homogeneous linear system are then the eigenvalues and eigenvectors of K, a set of 2n structural invariants of K. Importantly, λ_j can be complex as well as real.

In mechanics, modal analysis is typically concerned with vibrations and resonances in mechanical structures (Leonid Manevitch 2005). When their models are linear, or NL and linearized about their equilibrium solutions – a common practice in analyzing stability properties – some λ_j are complex. This means there exist oscillatory components in their dynamical behavior. In this context, linear modes are called linear normal modes[5], or LNMs, and the governing equations of motion can be decoupled into LNMs, e.g., by diagonalization of their system matrix. In decoupled form, i.e., with basis chosen to decouple the state variables, a linear system vibrates as if it consisted of independent oscillators governed by the (complex) modes, a property not shared with NL systems. Linear modal solutions also obey the superposition principle, i.e., free and forced oscillations are conveniently expressed as linear combinations of individual LNM motions in phase space. Importantly, modal solutions have an invariance property, i.e., a motion initiated on one specific LNM – representable as a two-dimensional planar surface in phase space – has no effect on the remaining LNMs, because they’re decoupled. This property permits the modal concept to be readily extended to NL systems.

By analogy, for NL models of the form

\dot{x} = f (x, u)

, a nonlinear normal mode (NNM) is defined similarly, as a two-dimensional invariant manifold in phase space. Nonlinear systems, however, can exhibit highly complex behaviors not possible with linear systems. These include (non-separable) modal interactions, saturation, limit cycles, bifurcations, jumps, subharmonic, superharmonic and internal resonances, and chaos. For this reason, NNM analysis requires more complex mathematical handling (Kerschen, Peeters et al. 2009).

Nonlinear modal analysis is about generating minimal reduced-order models that accurately capture the dynamics of higher-order models, with the help of NNMs. An effective and computationally efficient methodology for accomplishing this is based on the two-dimensional invariant-manifold approach, as developed in (Pierre, Jiang et al. 2006), for example. Geometrically, (NL) NNMs – like LNMs – are represented by two-dimensional nonplanar surfaces in phase space. They are, however, tangent to the planar surfaces of their linearized ODE LNMs, at the equilibrium point of their NL ODE system.

For mechanical system applications, the nonplanar surface NNM is parameterized by a single pair of state variables (displacement and velocity), chosen as master coordinates, with the remaining variables being functionally related to the chosen pair. In this special form, the original system behaves on the manifold surface like a nonlinear single-degree-of-freedom system. Solutions beginning in the manifold remain in it for all time, which extends the invariance property of LNMs to nonlinear systems[6].

Nonlinear Modes in Systems Biology

Conceptually, this treatment of NL normal modes in mechanics, as nonplanar two-dimensional manifolds in phase space, is directly extensible to NL systems biology models. The analogous application is reducing complex biosystem models using conservation relations, combined with dynamical approximations common in biochemical reactions. For example, enzymatic and cellular protein interaction networks (PIN) typically involve complex interactions among numerous cellular components and their models typically involve sets of coupled first-order ODEs, each representing the state variable motions of the interacting cellular components. Dominant characteristics of these systems is the presence of numerous inter-variable conservation relations, which serve to reduce system order, as well as numerous higher-order interactions that obey the quasi-steady state approximation (QSSA) (Segel and Slemrod 1989) or total quasi-steady state (tQSSA) assumptions (Tzafriri and Edelman 2004, Ciliberto 2007). In each of these latter interactions, the dominant dynamics are second-order, one each for the substrate and product, with all enzymatic reactions or other protein or metabolic interactions substituted by algebraic equations, based on the QSSA or tQSSA (Distefano 2014). On a microscopic level, each interaction obeying one of these assumptions also has two nonlinear modes, a fast one – on the “fast manifold” – and a slow one – on the “slow manifold”. The fast one disappears nearly instantaneously (or very quickly) relative to the slow one – which represents the relatively longer-term behavior of the mechanisms that dominate overall biosystem dynamics. A single mode, the slow one, strongly dominates the dynamics and approximates the subsystem interaction reasonably well. For example, this behavior is illustrated in the phase space figures (with nullcline “manifolds”) associated with a 4-dimensional Slyke-Cullen reaction feedback system model, where all trajectories starting at different points in phase space rapidly converge to a planar surface (Flach 2006).

This is a general extension of the QSSAs to higher-order biological systems, and is directly analogous to model reduction of NNM analysis of invariant manifolds on the phase space of mechanical systems. This means NNM analysis is equivalent to QSSA analysis, which likely came first (Segel and Slemrod 1989, Shaw 1994), with invariant manifolds of the biosystem in the substrate-product phase space of each reaction submodel.

Discussion

We’ve presented a cohesive viewpoint of the modal concept, merging classical linear with modern nonlinear dynamic system geometries. We used the classical multicompartmental modeling paradigm for developing this viewpoint – with compartments representing state variables (nodes), material flux functions or constant parameters (edges) connecting them, and I-O ports explicitly stated, all in accord with hypothesized real system structure. This has been presented as a unifying paradigm for establishing minimal homeomorphic model structure from I-O data signatures – represented as modal responses on 2-dimensional invariant manifolds – for NL as well as linear systems. Naturally encompassing these notions, compartmental modeling theory has a sound physical and mathematical basis, much of it motivated early on by problems in chemical kinetics (Hevesey 1923, Atkins 1969, Atkins 2005). Ironically, the new systems biology has had similar motivational forces – from biochemistry – but about a decade ago, the new theory separated itself from well-developed compartmentalization concepts. Developers of the Systems Biology Markup Language (SBML), a language designed for developing universally compatible software in biomodeling, ignored the past – likely inadvertently – in an important software “venue” receiving much attention in the systems biology modeling community at the time[7]. They adopted an older and looser notion of compartment than that honed and refined by the established community of 20th-century compartmental analysts. Mones Berman, John Jacquez and Sol Rubinow, for example, wrote books, ran national or university programs, or developed fairly widely used modeling software (SAAM etc) based on or about compartmental analysis methods – with compartments defined quite differently than in SBML lingo.

In the lingo of compartmental analysis, reviewed earlier, a compartment is an amount of material X that acts kinetically in a homogeneously distinct manner – ultimately a state variable in the language of dynamic systems and differential equations. In contrast, compartment in the “lingua franca” of SBML means a space that (typically) contains multiple state variables[8]. This is régressif (pardon my French), taking us back to an earlier time in biology when physical spaces, like organs in a biological system, were conveniently called compartments. Colloquially speaking, they still are to some, but more likely they contain numerous sub-compartments in a myriad of forms – like nuclear, or cytosolic, or tubule X, for example. As a consequence, both free and commercial software packages for systems biology modeling (SimBiology, Copasi, Cell Designer, etc) are written or being written with compartments that speak SBML, rather than classical modeling science, and many users (like myself) need an attitude adjustment as well as a workaround to use them.

Change is often useful, leading to progress, but in my view not in this case. I’ve done my best to show here that there’s good reason to retain the classical notion of compartment in dynamic systems biology modeling. In a nutshell, compartment as state variable promotes the congruence of relevant mechanistic model topology (connectedness) with modal I-O data. Classical compartments circumscribe and embody the deep connection between signal data (measurements) and model structure – linear or nonlinear – at the same time leveraging well-developed tools of engineering and physics – and compartmental modeling theory – based on ordinary differential equations.

Appendix – Pseudo Code for Calculating Mode Number for each Measurable Compartment[9]

Require:Graph G(V,E), node (compartment) i receives input

Ensure: Modes[j] mode-count for all compartments, for all j ϵ V

// Determine compartments reachable from the input

whilev ϵ V

if v is reachable from i then// Using DFS or BFS:

TrueMode[v] ← 1

end if

end while

// Count modes based on connectivity:

while v ϵ V do

while v’ ϵ V do

if v is reachable from v’ then // Using DFS or BFS:

if TrueMode[v] == 1 then

Modes[v] ← Modes[v] + 1

end if

end while

This algorithm actually finds the maximum number of modes in each compartment, because it is completely symbolic – with no regard for numerical values of model parameters. Hidden modes can be missed, e.g., those due to model parameters coinciding, like the k₂₂ = k₃₃ in Example 4.

Notes

1	Hidden modes (state variables, compartments) are directly analogous to hidden states in a Markov model (hidden Markov model).
2	All well and good. There would be no nomenclature incompatibility problem (new vs. old) if the differences stopped there, because species = state variable, or classical compartment = state variable, would amount to simply renaming things. Chemical species and compartment appear to be one and the same. But the SBML documention, amplifying on the definition, adds that species refers to ‘a pool of entities considered indistinguishable from each other for the purposes of the model, participate in reactions, and are located in a specific compartment (classical ‘space’ or container). Their pool then is our classical (homogeneous) compartment; and their (differently defined) ‘compartment’ is a ‘bounded space in which species are located’ (i.e. container of state variables). As noted above, in classical terms, pool is a superset of compartment, meant to represent an entity possibly distributed inhomogeneously in a distinct space (e.g. a biomolecular species in a poorly mixed container). They’ve redefined pool as well as compartment, with potential for confusion. The definitions of biostructural components chosen by molecular systems biologists are thus irreconcilable with the classical definitions in current use in other areas of life sciences. The meanings of species, compartment and pool-taken together-generates potential ambiguities for modeling.
3	This algorithm was originally written by Farhad Hormozdiari, and augmented by Teaching Assistant Long Nguyen, both PhD students at UCLA at the time of this writing.
4	Multicompartmental models also can have hidden oscillations (hidden oscillatory modes) in some compartments – with complex eigenvalues (modes) hidden from outputs measured in a central compartment. This interesting phenomenon can occur in potentially physiologically realizable systems, modeled with generalized mammillary models, as developed in (Fagarasan and DiStefano III 1986).
5	https://en.wikipedia.org/wiki/Normal_mode
6	More general theory and applications of NNM analysis are available in (Kerschen, Peeters et al. 2009).
7	http://en.wikipedia.org/wiki/SBML
8	Defined on the same site: http://en.wikipedia.org/wiki/SBML or the source site: http://sbml.org
9	This algorithm was written by Farhad Hormozdiari and Long Nguyen, PhD students at UCLA at the time of this writing. (will be made an Acknowledgement)

References

Alon, U. (2007). An Introduction to Systems Biology. Boca Raton, FL, Chapman & Hall/CRC.
Atkins, G. L. (1969). Multicompartment Models for Biological Systems. London, Methuen.
Atkins, P. , Jones, L. (2005). Chemical Principles: The Quest for Insight. New York, W. H. Freeman and Company.
Bergner, P. (1967). The concepts of mass, volume, and concentration. Compartments, Pools, and Spaces in Medical Physiology. P. Bergner and C. Lushbaugh. Washington, DC, US Atomic Energy Commission, Division of Technical Information: 21-37.
Chen, C. (1970). Introduction to Linear System Theory. New York, Holt, Rinehart, & Winston.
Chen, C. (1985). Introduction to Linear System Theory. New York, Holt, Rinehart, & Winston.
Ciliberto, A. , Capuani, F., Tyson, J. (2007). "Modeling Networks of Coupled Enzymatic Reactions Using the Total Quasi-Steady State Approximation." PLoS Comput Biol 3(3): e45.
Cormen, T. H. , Leiserson, C., Rivest,R. L. and Stein, C. (2009). Introduction to Algorithms. Cambridge, MIT Press.
DiStefano III, J. J. (2014). Dynamic systems biology modeling and simulation, Academic Press/Elsevier.
Distefano, J. (2019). Dynamic Biosystem Modeling & Simulation Methodology: Integrated & Accessible, Amazon-KDP.
DiStefano, J. J. I. (2013). Dynamic systems biology modeling and simulation, Academic Press/Elsevier.
Fagarasan, J. and J. DiStefano III (1986). "Hidden pools, hidden modes, and visible repeated eigenvalues in compartmental models." Mathematical Biosciences 82: 87-113.
Flach, E. H. , Schnell, S. (2006). "Use and abuse of the quasi-steady-state approximation." IEE Proc.-Syst. Biol. 153(4): 187-191.
Godfrey, K. (1983). Compartmental Models and their Application. New York, Academic Press.
Hevesey, G. (1923). "The absorption and translocation of lead by plants.
J. Applied Physiol. 17: 439-445.
Jacquez, J. A. (1996). Compartmental Analysis in Biology and Medicine. Ann Arbor, Biomedware.
Kerschen, G., M. Peeters, J. C. Golinval and A. F. Vakakis (2009). "Nonlinear normal modes, Part I: A useful framework for the structural dynamicist." Mechanical Systems and Signal Processing 23(1): 170-194.
Klipp, E. , Liebermeister, Wolfram, Wierling, Christoph, Kowald, Axel, Lehrach, Hans,Herwig, Ralf (2009). Systems Biology: A Textbook, Wiley-VCH, Weinheim.
Leonid Manevitch, A. I. M. (2005). The Mechanics Of Nonlinear Systems With Internal Resonances, Imperial College Press.
Palsson, B. (2006). Systems Biology: Properties of Reconstructed Networks. Cambridge, Cambridge University Press.
Pierre, C., D. Jiang and S. Shaw (2006). "Nonlinear normal modes and their application in structural dynamics." Mathematical Problems in Engineering 2006.
Rescigno, A. , Thakur, A., Bertrand,K, Brill, A, Mariani, G (1990). "Tracer kinetics: a proposal for unified symbols and nomenclature." Phys. Med. Biol. 35(3): 449-465.
Rubinow, S. (1975). Introduction to Mathematical Biology. New York, John Wiley.
Rubinow, S. (1975). "On closed or almost closed compartment systems." Mathematical Biosciences 18: 245-253.
Segel, L. A. and M. Slemrod (1989). "The Quasi-Steady-State Assumption: A Case Study in Perturbation." SIAM Review 31(3): 446-477.
Shaw, S. , Pierre, C. (1994). "Normal Modes of Vibration for Non-Linear Continuous Systems." J. Sound and Vibration 169(3): 319-347.
Tzafriri, A. R. and E. R. Edelman (2004). "The total quasi-steady-state approximation is valid for reversible enzyme kinetics." Journal of Theoretical Biology 226(3): 303-313.
Voit, E. (2012). A First Course in Systems Biology, Garland Science.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer