# Effect of Linearization in a WNT Signaling Model

Authors: C. F. Ciușdel; S. Coman; Cr. Boldișor; T. Kessler; A. Muradyan; A. Kovachev; H. Lehrach; C. Wierling; L. M. Itu
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

A nonlinear model consisting of a system of coupled ordinary differential equations (ODE), describing a biological process linked with cancer development, is linearized using Taylor series and tested against different magnitudes of input perturbations, in order to investigate the extent to which the linearization is accurate. The canonical wingless/integrated (WNT) signaling pathway is considered. The linearization procedure is described, and special considerations for linearization validity are analyzed. The analytical properties of nonlinear and linearized systems are studied, including aspects such as existence of steady state and initial value sensitivity. Linearization is a useful tool for speeding up drug response computations or for providing analytical answers to problems such as required drug concentrations. A Monte Carlo-based error testing workflow is employed to study the errors introduced by the linearization for different input conditions and parameter vectors. The deviations between the nonlinear and the linearized system were found to increase in a polynomial fashion w.r.t. the magnitude of tested perturbations. The linearized system closely followed the original one for perturbations of magnitude within 10% of the base input vector which yielded the state-space fixed point used for the linearization.


---

## Body

## 1. Introduction

The WNT signaling pathway is a cascade of biochemical reactions transducing signals from the extracellular space into the cells. This pathway is involved in different biological processes such as embryonic development, tissue homeostasis, and tumorigenesis [1–3]. The pathway is activated by the binding of WNT-protein ligands to Frizzled family receptors, which pass the biological signal to the Dishevelled (DVL) protein inside the cell.Specifically, we considered the canonical WNT pathway that regulates translocation of the cytoplasmicβ-catenin (CTNNB1) into the nucleus. Without WNT signaling, β-catenin is degraded by a destruction complex, which includes AXIN, APC, GSK3β, and CK. The last two components phosphorylate β-catenin on several serine and threonine residues targeting it for ubiquitination and subsequent proteosomal degradation. This continual elimination of the β-catenin prevents its translocation to the nucleus [1–3].The pathway is activated when WNT ligands bind FZD and LRP5/6. This disrupts the function of the destruction complex by recruiting it to the plasma membrane, leading toβ-catenin stabilization in the cytoplasm and its subsequent nuclear translocation. In the nucleus, β-catenin binds TCF/LEF transcription factor proteins and activates expression of WNT target genes [3].Nuclearβ-catenin together with 5 different TCF/LEF transcription factors and 4 other proteins (BCL9, BCL9L, PYGO1, and PYGO2) makes 20 positive readouts. Without WNT signaling and, corresponding nuclear β-catenin, the 5 TCF/LEF transcription factors are combined with 4 transcriptional corepressor proteins forming 20 negative readouts, thus leading to a total of 40 readouts. Figure 1 depicts an overview on the WNT model structure.Figure 1
Overview of the WNT signaling model. The WNT signaling pathway model contains 7 WNT ligands, 7 FZD receptors, 2 LRP coreceptors, and 40 readouts.The main objective of this paper is to investigate the effect of linearization (based on Taylor series expansion) for a custom prototype model of the human WNT signaling pathway described above, similar to the Reactome model in terms of structure and complexity [4, 5]. The approach considered herein is however generic and applicable to any nonlinear system of coupled ordinary differential equations (ODE).Linearization may be a useful approach for many signaling model-based use cases. First of all, we refer to the case where only “relative” data are available. “Relative” datasets do not contain pairs of input-output vectors of absolute values. Instead, a data sample represents the ratio of outputs for two different input conditions/vectors. The datasets may be organized in chunks. A chunk has a unique base input vector; other input vectors are variations of the base one. Each input vector has a corresponding output ratio, i.e., the ratio between its corresponding output vector and the base output vector. These ratios may be more relevant than plain absolute values, as they can easily indicate up or down regulation of model components/states for various input conditions/vectors.To evaluate a ratio, the steady state corresponding to the base input vector must be computed first. Using the base steady-state point as initial conditions, other forward simulations compute the steady state for each input vector. This procedure usually employs numeric integrators for ODE systems, which have various computational and time requirements. A linearization may be conducted around the base steady-state point. For inputs which do not lead to large linearization errors (i.e., ones sufficiently close to the base input vector), the linearized model can be employed since it is faster and easier to use than the nonlinear one.To compute the numerical solution of stiff nonlinear ODE systems (until steady state is reached) a system of linear equations needs to be solved multiple times. Depending on the time step, and on the number of time steps required to reach the steady state, the total number of solutions may vary between a few tens and several thousands. In contrast, for the linearized system, a system of linear equations must be solved exactly once. Furthermore, additional runtime is required for evaluating the model equations and the Jacobian.Another advantageous use case for linearization is the computation of the required input deviations (relative to the base input vector) that would lead to a desired change in the state variable vector (relative to the base steady-state point). This approach may be employed to compute the ideal specificity of a potential drug or drug combination to treat a certain type of tumor, or to help in designing new drugs for specific patients or groups of patients.Moreover, linearization is a step often performed [6–8] in designing control systems for regulating different components of a biological system, either through a nonlinear feedback, or by linearizing the model’s equations, as it allows for a direct application of both classical and modern control strategies.Studying the accuracy of the linearization is therefore useful for indicating when the linearized model can be successfully employed, i.e., when the linearization errors are negligible. These errors depend on the distance between the base input vector and other input vectors. If a suitable threshold can be found, input vectors meeting the threshold can be analyzed using the linearized model, while the rest can be forward simulated using the nonlinear model.In the following, the base input vector is denoted asub, and the other input vectors are regarded as input perturbations.


## 2. Methods

### 2.1. Nonlinear Systems and Linearization

Generically, a dynamic nonlinear system is described using the following equations:(1)dxdt=fx,p,u,(2)y=gx,p,u,where x is the state vector, p is the parameter vector, u is the input vector, and y is the output vector. Function f describes the behavior of the system, while function g defines the output variables of the model. If f is Lipschitz continuous, then equation (1) has a unique solution [9].For the considered WNT model,f and g are Lipschitz continuous functions and parameter vector p is presumed constant (therefore, its notation can be omitted). The model has the following size: x∈ℝ421, u∈ℝ60, p∈ℝ1154, and y∈ℝ40.If input vectoru has a step-like evolution, then at a fixed point xss (considered as steady state), equations (1) and (2) are rewritten as follows:(3)0=fxss,ubf,(4)yss=gxss,ubf,where yss is the steady-state output vector and ubf is the final value of the base input vector ub.If a linearization is performed aroundxss, using Taylor series expansion, and all terms of second or higher order are being discarded, then(5)flinearx,u=fxss,ubf+∂f∂xxss,ubfx−xss+∂f∂uxss,ubfu−ubf,where ∂f/∂xxss,ubf is the state Jacobian evaluated at xss,ubf and ∂f/∂uxss,ubf is the input Jacobian evaluated at xss,ubf. Vectors x−xss and u−ubf are states and input vector deviations from the point where the linearization was performed at.Typically, the following notations are used in control theory [10, 11]:(6)A=∂f∂xxss,ubf,(7)B=∂f∂uxss,ubf.If conditions (i-ii) listed below hold, the nonlinear system is already linear in its input space and the linearization does not introduce any errors from omitting nonlinear terms ofu−ubf since they do not exist. These conditions are true for the WNT model considered in this paper:(i)
The input Jacobian is constant (i.e., does not depend onx or u)(ii)
The state Jacobian does not depend explicitly onuUsing equations (3) and (6)-(7),(8)dxlineardt=flinearxlinear,u=Axlinear−xss+Bu−ubf,where xlinear is the state vector of the linearized system.Based on equation (3), at the point of linearization, dxss/dt is zero. Therefore,(9)ddtxlinear−xss=Axlinear−xss+Bu−ubf.Ifxdev is defined as xlinear−xss, i.e., the vector representing the deviations of the state from the point of linearization for the linear model, then (10)dxdevdt=Axdev+Bu−ubf.The agreement between the linearized model and the original model can be formalized as follows:(11)dxdevdt=Axdev+Bu−ubf,dxnonlindt=fxnonlin,u,e=xnonlin−xdev−xss,where dxnonlin/dt is the state evolution as dictated by the original nonlinear system and e is the error vector representing the differences introduced by the linearization. A new nonlinear system was obtained, describing the evolution in time of the agreement between the linearized and the original model. By defining(12)xagreement=xdevxnonlin,an analysis related to the errors introduced by the linearization can be conducted for different input vectors u. An example workflow would be as follows:(a)
Forward simulate the nonlinear model until a steady-state pointxss is reached for input ub(b)
Compute matricesA and B(c)
Forward simulate model (11) for a new input vector u, with the initial conditions xagreement0=0xss(d)
Compute measures of interest on the agreement error vectore


### 2.2. Effects of Linearization

If the linearization is performed around a hyperbolic fixed point, the Hartman–Grobman theorem [12] guarantees that the linearized system will exhibit the same qualitative behavior as the original system, in a sufficiently small neighborhood around the linearization point. A fixed point xf is “hyperbolic” if all eigenvalues of the Jacobian evaluated at xf (i.e., the linearized systems’ poles) have a nonzero real part (i.e., the system does not have purely oscillatory modes).In terms of quantitative behavior, usually a linearization is considered accurate only for small deviations of the state vector around an equilibrium point. Of importance is the steady-state behavior of the linearized system when compared to the original one, i.e., the steady-state deviations/errors (introduced by the simplified model) for different magnitudes of exogenous input perturbations.In case of stable linearized models, for step-like inputs, the steady state is dependent only on the final value of the exogenous input vector and on the model matrices:(13)0=dxdevssdt=Axdevss+But⟶∞−ubf,xdevss=−A−1But⟶∞−ubf.Herein, all vector fieldsu which have the following properties are regarded as “step-like” inputs:(14)ut⟶∞=limt⟶∞ut=ufinal=const,ufinal≠ubf.This translates to the new steady state depending neither on the time-dependent evolution of the input vector (it depends only on its final value) nor on the initial conditionsx0.However, these considerations are not universally true for nonlinear systems. The existence of the steady state is guaranteed if the following lemma applies. SupposeB is a nonempty bounded subset of ℝn. By definition [13], the ω limit of set B, denoted ωB, is the set of all points x for which a sequence of pair xk,tk exists, with xk∈B and limk⟶∞tk=∞ such that(15)limk⟶∞xtk,xk=x.Lemma.
[13]. Let B be a nonempty bounded subset of ℝn and suppose there is a number M such that xt,x0≤M for all t≥0 and all x0∈B. Then, ωB is a nonempty compact invariant set. Moreover, the distance of xt,x0 from ωB approaches 0 as t⟶∞, uniformly in x0∈B.Less formally stated, if the nonlinear system has an overall stable behavior, i.e., with initial conditionsx0 in a neighborhood B, its state does not diverge towards infinity, then at a large enough time t≥T, its state trajectory will enter and remain in a closed set ωB. If neighborhood B is centered around an attractor fixed point xss, which can be used as a linearization point, the lemma suggests, if its conditions are met, that both the linearized and the nonlinear system will reach a new steady state in case of an exogenous perturbation. It is worth noting that the lemma deals explicitly with autonomous systems, but, if the evolution of the exogenous perturbation is known, it can be integrated into an extended autonomous system, comprising two subsystems: the perturbation model and the original nonlinear model [14]. In the following, let ωB consist of only one point, i.e., the system has no oscillatory behavior at steady state.To analyze the influence of the initial conditionx0 on the state trajectory of the nonlinear model, the Lyapunov exponents theory can be employed. We consider Φ as the discrete-time counterpart of the continuous-time model f, with the property that xk+1=Φxk, where xk=xt0+kTs and Ts is the sampling period. The following statement holds (any conclusions drawn for the discrete model are also valid for the continuous one): the global Lyapunov exponent of Φ at x with respect to direction y is defined as the limit (if it exists) [15]:(16)limt⟶∞1tlimδ⟶0lnΦtx+δy−Φtxδy=λ,where Φkx0=xtk,x0. The global Lyapunov exponent λ describes the decay in time of the distance between two trajectories which start at two separate initial points x and x+δy:(17)Φtx+δy−Φtxδy≈eλt.Ifλ<0, any state trajectory which starts at x+δy (for any sufficiently small δ) will converge to the “main” trajectory which starts at x. Since for a small δ we can state that(18)Φx+δy−Φx≈JΦ,xδy,where JΦ,x is the Jacobian of Φ evaluated at x and λ can be expressed as(19)λ=limn⟶∞1nlnJΦn,xyy=1n∑t=0n−1lnJΦt+1,xyJΦt,xy=1n∑t=0n−1lnJΦ,ΦtJΦt,xyJΦt,xy.Thus, the global Lyapunov exponent ofΦ at point x with respect to direction y is the average of the local Lyapunov exponents of Φ at points x0,x1,x2,… of the trajectory of x w.r.t direction JΦk,xy at each point xk [15]. The local Lyapunov exponents are characterized by the eigenvalues of JΦ,xk. Therefore, if the system reaches a new steady state (e.g., an attractor fixed point xssfinal) when an input perturbation is applied (starting from an initial stable equilibrium point xss0 which can be employed for the linearization), all state trajectories starting from a sufficiently small neighborhood B of xss0 will converge to xssfinal.LetS⊆ℝn be a compact and convex set of points x for which the state Jacobian evaluated at x has negative real-part eigenvalues. More formally, let λi be the eigenvalues of ∂f/∂xx,p, and ℜλi<0,∀x∈S,i∈1,…,n. In this case, the eigenvalues of the discrete-time model Φ lie inside the unit circle. S is restricted so that Φx∈S,∀x∈S. Then, at each point x∈S, the local Lyapunov exponent w.r.t. direction y is defined as(20)limδ⟶0lnΦx+δy−Φxδy,and is negative; therefore,(21)Φx+Δy−Φx≤cΔy,where 0<c<1 and Δ>0 is a small enough number.It can be shown that, for a constant inputu, for every initial point x0∈S, the nonlinear model reaches a unique fixed point xss∈S at steady state. Banach’s fixed point theorem states that if Φ is a strict contraction, then Φ has a unique fixed point in S [9]. A map Φ is a strict contraction on S, if ∀p,q∈S:(22)Φp−Φq≤cp−q,where c<1. To prove equation (22), consider the line segment connecting points p and q, and the sequence of distinct points p1,p2,…pk,qk,qk−1,…,q1 between p and q. Then,(23)Φp−Φq=Φp−Φp1−Φq−Φq1+Φp1−Φq1≤≤Φp−Φp1+Φq−Φq1+Φp1−Φq1.The points of each pair (pk−1, pk) are close enough:(24)Φpk−1−Φpk≤cpkpk−1−pk,where cpk<1 is the local Lyapunov exponent of Φ at pk. Applying the procedure of equation (23) iteratively and employing equation (24) yields(25)Φp−Φq≤cp1p−p1+…+cpk−1pk−pk−1+…+cq1q−q1.Letcmax=maxcp1,…,cpk−1,…,cq1<1.Then,(26)Φp−Φq≤cmaxp−p1+…+pk−pk−1+…+q−q1.Becausep,p1,p2,…pk,qk,qk−1,…,q1,q are collinear points, their inner distances sum to p−q. Therefore,(27)Φp−Φq≤cmaxp−q,and Φ is a strict contraction for all p,q∈S. This also implies that the ω limit of S consists of only one point. For comparison, a stable linear system will always reach a unique steady-state fixed point (i.e., ∀x0∈ℝn at t0) for a specific constant input, as stated in equation (13).As the global Lyapunov exponent is the average of the local Lyapunov exponents, the existence ofS is a particular case of the global Lyapunov exponent reasoning described above. S may also be regarded as a smaller subset of the attraction basin of the attractor fixed point xss.In conclusion, a linearization offers valid results only when tested against valid input perturbations, i.e., perturbations that do not drive the nonlinear system into an unstable region. In this case, the nonlinear and linearized models have similar behaviors: existence of steady state and (local) insensitivity w.r.t. initial conditions. The application of an exogenous input perturbation modifies the location in the state space of the attractor fixed pointxss, and therefore, it invalidates the prerequisites of the Hartman–Grobman theorem (due to the fact that the linearized system was defined at the initial fixed point location xss0, before applying the exogenous perturbation). Although the trajectories of the models will initially behave qualitatively the same when starting from xss0, i.e., the trajectories will start varying in the same direction (since the linearization is accurate for small deviations around xss0), the two models may behave qualitatively differently around xss. This holds true especially if the Jacobian varies considerably between xss and xss0.The actual steady-state error between the linearized model and the original model is influenced by the Hessian of the system. For exemplification, suppose conditions (i-ii) hold, letxlin be an attractor fixed point around which a linearization was performed, and the input perturbation u be a Heaviside step function. As a first example, consider the Taylor expansion of f around xlin:(28)fx,u=fxlin,ubf+∂fxlin∂xx−xlin+12!x−xlinT∂2fxlin∂x2x−xlin+…+Bu−ubf.At a new steady state of the nonlinear modelxss, fxss,u will be zero. If the Hessian norm is relatively large, the large weights of the second-order term will introduce errors between the steady state of the linearized and the steady state of the original model.As a second example, computing the partial derivatives ofx w.r.t. u at both fixed points xlin and xss yields(29)∂xlin∂u=−∂fxlin,u∂x−1B=−A−1B,∂xss∂u=−∂fxss,u∂x−1B.If the Jacobian varies significantly between the two fixed pointsxlin and xss, the linearized and the original system will exhibit different input sensitivities. The rate of change of the Jacobian w.r.t. δx≈xss−xlin is dictated by the Hessian.As a third example, by following an iterative procedure, one can forward simulate the linearized system for small enough time increments so that it follows precisely the nonlinear system. More formally, there is aδ>0 so that, for every ε>0, if xnonlineart−xlineart<ε, then Δt=tk+1−tk<δ. At the end of each time increment, the linearization (i.e., the Jacobian) is recomputed. Applying this for k steps until the new steady state is reached yields(30)0=fxss,u=∑k∂fxk∂xxk+1−xk+Bu−ubf.If the Jacobian is constant, equation (30) reduces to equation (13). Therefore, if the Jacobian varies significantly along the state trajectory towards the new attractor point given by the applied perturbation, the linearized system will deviate from the real one. Again, if the Hessian computed at xlin has a relatively large norm, linearization errors are to be expected.


### 2.3. Properties of the Considered WNT Model

By analyzing the model equations, we observe that functionf has for all state variables the following form:(31)x˙k=fkx,p,u=−xkh1kxi,p+h2kxi,p+uk,where xk is the kth state variable, xi is the set of all state variables excluding xk (i.e., i∈1,…,k−1,k+1,…,421), and h1k and h2k are nonlinear Lipschitz continuous functions with the property that h1kx,p>ck>0 and h2kx,p≥0, ∀x,p>0, and ck is constant. When both the exogenous input uk and the parameter vector p are chosen to be positive (which is the case for the considered datasets) and all state variables are nonnegative, the subsystem defined by fk is equivalent to a stable linear time-variant first-order filter. This, however, does not imply that the entire model is stable, due to the existence of nonlinear feedback loops.If all state variables have nonnegative initial values (att0), xk will always remain nonnegative. To prove this, let all other state variables be nonnegative and consider the zerocrossing of xk. Equation (31) is reduced to(32)x˙k=fkx,p,u=h2kxi,p+uk≥0.Therefore,xk will increase again as a positive number. As a consequence, negative values have no relevance for the nonlinear model and can simply be corrected in the linearized model. The predictions of the linearized model were thus in the end adjusted by clipping the state variable values in the range 0,∞.


### 2.4. Monte Carlo Analysis of Linearization Accuracy Range

To estimate the quantitative behavior of the linearized model, the following workflow was used:(1)
Generate a random parameter vectorp and set rnoise, representing the range of a uniform distribution centered around 1, used to perturb the input vectors(2)
For each input sampleub,i,(a)
Forward simulate the nonlinear system until the steady-statexss0 is reached(b)
Linearize the nonlinear model aroundxss0(c)
Computeui by perturbing input vector ubf,i, with a multiplicative noise of range ±rnoise%(d)
Forward simulate system (11) until the steady state, using the initial conditions xagreement0=0xss0(3)
Compute metrics by averaging over all input samples, for the currentp and rnoise.


## 2.1. Nonlinear Systems and Linearization

Generically, a dynamic nonlinear system is described using the following equations:(1)dxdt=fx,p,u,(2)y=gx,p,u,where x is the state vector, p is the parameter vector, u is the input vector, and y is the output vector. Function f describes the behavior of the system, while function g defines the output variables of the model. If f is Lipschitz continuous, then equation (1) has a unique solution [9].For the considered WNT model,f and g are Lipschitz continuous functions and parameter vector p is presumed constant (therefore, its notation can be omitted). The model has the following size: x∈ℝ421, u∈ℝ60, p∈ℝ1154, and y∈ℝ40.If input vectoru has a step-like evolution, then at a fixed point xss (considered as steady state), equations (1) and (2) are rewritten as follows:(3)0=fxss,ubf,(4)yss=gxss,ubf,where yss is the steady-state output vector and ubf is the final value of the base input vector ub.If a linearization is performed aroundxss, using Taylor series expansion, and all terms of second or higher order are being discarded, then(5)flinearx,u=fxss,ubf+∂f∂xxss,ubfx−xss+∂f∂uxss,ubfu−ubf,where ∂f/∂xxss,ubf is the state Jacobian evaluated at xss,ubf and ∂f/∂uxss,ubf is the input Jacobian evaluated at xss,ubf. Vectors x−xss and u−ubf are states and input vector deviations from the point where the linearization was performed at.Typically, the following notations are used in control theory [10, 11]:(6)A=∂f∂xxss,ubf,(7)B=∂f∂uxss,ubf.If conditions (i-ii) listed below hold, the nonlinear system is already linear in its input space and the linearization does not introduce any errors from omitting nonlinear terms ofu−ubf since they do not exist. These conditions are true for the WNT model considered in this paper:(i)
The input Jacobian is constant (i.e., does not depend onx or u)(ii)
The state Jacobian does not depend explicitly onuUsing equations (3) and (6)-(7),(8)dxlineardt=flinearxlinear,u=Axlinear−xss+Bu−ubf,where xlinear is the state vector of the linearized system.Based on equation (3), at the point of linearization, dxss/dt is zero. Therefore,(9)ddtxlinear−xss=Axlinear−xss+Bu−ubf.Ifxdev is defined as xlinear−xss, i.e., the vector representing the deviations of the state from the point of linearization for the linear model, then (10)dxdevdt=Axdev+Bu−ubf.The agreement between the linearized model and the original model can be formalized as follows:(11)dxdevdt=Axdev+Bu−ubf,dxnonlindt=fxnonlin,u,e=xnonlin−xdev−xss,where dxnonlin/dt is the state evolution as dictated by the original nonlinear system and e is the error vector representing the differences introduced by the linearization. A new nonlinear system was obtained, describing the evolution in time of the agreement between the linearized and the original model. By defining(12)xagreement=xdevxnonlin,an analysis related to the errors introduced by the linearization can be conducted for different input vectors u. An example workflow would be as follows:(a)
Forward simulate the nonlinear model until a steady-state pointxss is reached for input ub(b)
Compute matricesA and B(c)
Forward simulate model (11) for a new input vector u, with the initial conditions xagreement0=0xss(d)
Compute measures of interest on the agreement error vectore


## 2.2. Effects of Linearization

If the linearization is performed around a hyperbolic fixed point, the Hartman–Grobman theorem [12] guarantees that the linearized system will exhibit the same qualitative behavior as the original system, in a sufficiently small neighborhood around the linearization point. A fixed point xf is “hyperbolic” if all eigenvalues of the Jacobian evaluated at xf (i.e., the linearized systems’ poles) have a nonzero real part (i.e., the system does not have purely oscillatory modes).In terms of quantitative behavior, usually a linearization is considered accurate only for small deviations of the state vector around an equilibrium point. Of importance is the steady-state behavior of the linearized system when compared to the original one, i.e., the steady-state deviations/errors (introduced by the simplified model) for different magnitudes of exogenous input perturbations.In case of stable linearized models, for step-like inputs, the steady state is dependent only on the final value of the exogenous input vector and on the model matrices:(13)0=dxdevssdt=Axdevss+But⟶∞−ubf,xdevss=−A−1But⟶∞−ubf.Herein, all vector fieldsu which have the following properties are regarded as “step-like” inputs:(14)ut⟶∞=limt⟶∞ut=ufinal=const,ufinal≠ubf.This translates to the new steady state depending neither on the time-dependent evolution of the input vector (it depends only on its final value) nor on the initial conditionsx0.However, these considerations are not universally true for nonlinear systems. The existence of the steady state is guaranteed if the following lemma applies. SupposeB is a nonempty bounded subset of ℝn. By definition [13], the ω limit of set B, denoted ωB, is the set of all points x for which a sequence of pair xk,tk exists, with xk∈B and limk⟶∞tk=∞ such that(15)limk⟶∞xtk,xk=x.Lemma.
[13]. Let B be a nonempty bounded subset of ℝn and suppose there is a number M such that xt,x0≤M for all t≥0 and all x0∈B. Then, ωB is a nonempty compact invariant set. Moreover, the distance of xt,x0 from ωB approaches 0 as t⟶∞, uniformly in x0∈B.Less formally stated, if the nonlinear system has an overall stable behavior, i.e., with initial conditionsx0 in a neighborhood B, its state does not diverge towards infinity, then at a large enough time t≥T, its state trajectory will enter and remain in a closed set ωB. If neighborhood B is centered around an attractor fixed point xss, which can be used as a linearization point, the lemma suggests, if its conditions are met, that both the linearized and the nonlinear system will reach a new steady state in case of an exogenous perturbation. It is worth noting that the lemma deals explicitly with autonomous systems, but, if the evolution of the exogenous perturbation is known, it can be integrated into an extended autonomous system, comprising two subsystems: the perturbation model and the original nonlinear model [14]. In the following, let ωB consist of only one point, i.e., the system has no oscillatory behavior at steady state.To analyze the influence of the initial conditionx0 on the state trajectory of the nonlinear model, the Lyapunov exponents theory can be employed. We consider Φ as the discrete-time counterpart of the continuous-time model f, with the property that xk+1=Φxk, where xk=xt0+kTs and Ts is the sampling period. The following statement holds (any conclusions drawn for the discrete model are also valid for the continuous one): the global Lyapunov exponent of Φ at x with respect to direction y is defined as the limit (if it exists) [15]:(16)limt⟶∞1tlimδ⟶0lnΦtx+δy−Φtxδy=λ,where Φkx0=xtk,x0. The global Lyapunov exponent λ describes the decay in time of the distance between two trajectories which start at two separate initial points x and x+δy:(17)Φtx+δy−Φtxδy≈eλt.Ifλ<0, any state trajectory which starts at x+δy (for any sufficiently small δ) will converge to the “main” trajectory which starts at x. Since for a small δ we can state that(18)Φx+δy−Φx≈JΦ,xδy,where JΦ,x is the Jacobian of Φ evaluated at x and λ can be expressed as(19)λ=limn⟶∞1nlnJΦn,xyy=1n∑t=0n−1lnJΦt+1,xyJΦt,xy=1n∑t=0n−1lnJΦ,ΦtJΦt,xyJΦt,xy.Thus, the global Lyapunov exponent ofΦ at point x with respect to direction y is the average of the local Lyapunov exponents of Φ at points x0,x1,x2,… of the trajectory of x w.r.t direction JΦk,xy at each point xk [15]. The local Lyapunov exponents are characterized by the eigenvalues of JΦ,xk. Therefore, if the system reaches a new steady state (e.g., an attractor fixed point xssfinal) when an input perturbation is applied (starting from an initial stable equilibrium point xss0 which can be employed for the linearization), all state trajectories starting from a sufficiently small neighborhood B of xss0 will converge to xssfinal.LetS⊆ℝn be a compact and convex set of points x for which the state Jacobian evaluated at x has negative real-part eigenvalues. More formally, let λi be the eigenvalues of ∂f/∂xx,p, and ℜλi<0,∀x∈S,i∈1,…,n. In this case, the eigenvalues of the discrete-time model Φ lie inside the unit circle. S is restricted so that Φx∈S,∀x∈S. Then, at each point x∈S, the local Lyapunov exponent w.r.t. direction y is defined as(20)limδ⟶0lnΦx+δy−Φxδy,and is negative; therefore,(21)Φx+Δy−Φx≤cΔy,where 0<c<1 and Δ>0 is a small enough number.It can be shown that, for a constant inputu, for every initial point x0∈S, the nonlinear model reaches a unique fixed point xss∈S at steady state. Banach’s fixed point theorem states that if Φ is a strict contraction, then Φ has a unique fixed point in S [9]. A map Φ is a strict contraction on S, if ∀p,q∈S:(22)Φp−Φq≤cp−q,where c<1. To prove equation (22), consider the line segment connecting points p and q, and the sequence of distinct points p1,p2,…pk,qk,qk−1,…,q1 between p and q. Then,(23)Φp−Φq=Φp−Φp1−Φq−Φq1+Φp1−Φq1≤≤Φp−Φp1+Φq−Φq1+Φp1−Φq1.The points of each pair (pk−1, pk) are close enough:(24)Φpk−1−Φpk≤cpkpk−1−pk,where cpk<1 is the local Lyapunov exponent of Φ at pk. Applying the procedure of equation (23) iteratively and employing equation (24) yields(25)Φp−Φq≤cp1p−p1+…+cpk−1pk−pk−1+…+cq1q−q1.Letcmax=maxcp1,…,cpk−1,…,cq1<1.Then,(26)Φp−Φq≤cmaxp−p1+…+pk−pk−1+…+q−q1.Becausep,p1,p2,…pk,qk,qk−1,…,q1,q are collinear points, their inner distances sum to p−q. Therefore,(27)Φp−Φq≤cmaxp−q,and Φ is a strict contraction for all p,q∈S. This also implies that the ω limit of S consists of only one point. For comparison, a stable linear system will always reach a unique steady-state fixed point (i.e., ∀x0∈ℝn at t0) for a specific constant input, as stated in equation (13).As the global Lyapunov exponent is the average of the local Lyapunov exponents, the existence ofS is a particular case of the global Lyapunov exponent reasoning described above. S may also be regarded as a smaller subset of the attraction basin of the attractor fixed point xss.In conclusion, a linearization offers valid results only when tested against valid input perturbations, i.e., perturbations that do not drive the nonlinear system into an unstable region. In this case, the nonlinear and linearized models have similar behaviors: existence of steady state and (local) insensitivity w.r.t. initial conditions. The application of an exogenous input perturbation modifies the location in the state space of the attractor fixed pointxss, and therefore, it invalidates the prerequisites of the Hartman–Grobman theorem (due to the fact that the linearized system was defined at the initial fixed point location xss0, before applying the exogenous perturbation). Although the trajectories of the models will initially behave qualitatively the same when starting from xss0, i.e., the trajectories will start varying in the same direction (since the linearization is accurate for small deviations around xss0), the two models may behave qualitatively differently around xss. This holds true especially if the Jacobian varies considerably between xss and xss0.The actual steady-state error between the linearized model and the original model is influenced by the Hessian of the system. For exemplification, suppose conditions (i-ii) hold, letxlin be an attractor fixed point around which a linearization was performed, and the input perturbation u be a Heaviside step function. As a first example, consider the Taylor expansion of f around xlin:(28)fx,u=fxlin,ubf+∂fxlin∂xx−xlin+12!x−xlinT∂2fxlin∂x2x−xlin+…+Bu−ubf.At a new steady state of the nonlinear modelxss, fxss,u will be zero. If the Hessian norm is relatively large, the large weights of the second-order term will introduce errors between the steady state of the linearized and the steady state of the original model.As a second example, computing the partial derivatives ofx w.r.t. u at both fixed points xlin and xss yields(29)∂xlin∂u=−∂fxlin,u∂x−1B=−A−1B,∂xss∂u=−∂fxss,u∂x−1B.If the Jacobian varies significantly between the two fixed pointsxlin and xss, the linearized and the original system will exhibit different input sensitivities. The rate of change of the Jacobian w.r.t. δx≈xss−xlin is dictated by the Hessian.As a third example, by following an iterative procedure, one can forward simulate the linearized system for small enough time increments so that it follows precisely the nonlinear system. More formally, there is aδ>0 so that, for every ε>0, if xnonlineart−xlineart<ε, then Δt=tk+1−tk<δ. At the end of each time increment, the linearization (i.e., the Jacobian) is recomputed. Applying this for k steps until the new steady state is reached yields(30)0=fxss,u=∑k∂fxk∂xxk+1−xk+Bu−ubf.If the Jacobian is constant, equation (30) reduces to equation (13). Therefore, if the Jacobian varies significantly along the state trajectory towards the new attractor point given by the applied perturbation, the linearized system will deviate from the real one. Again, if the Hessian computed at xlin has a relatively large norm, linearization errors are to be expected.


## 2.3. Properties of the Considered WNT Model

By analyzing the model equations, we observe that functionf has for all state variables the following form:(31)x˙k=fkx,p,u=−xkh1kxi,p+h2kxi,p+uk,where xk is the kth state variable, xi is the set of all state variables excluding xk (i.e., i∈1,…,k−1,k+1,…,421), and h1k and h2k are nonlinear Lipschitz continuous functions with the property that h1kx,p>ck>0 and h2kx,p≥0, ∀x,p>0, and ck is constant. When both the exogenous input uk and the parameter vector p are chosen to be positive (which is the case for the considered datasets) and all state variables are nonnegative, the subsystem defined by fk is equivalent to a stable linear time-variant first-order filter. This, however, does not imply that the entire model is stable, due to the existence of nonlinear feedback loops.If all state variables have nonnegative initial values (att0), xk will always remain nonnegative. To prove this, let all other state variables be nonnegative and consider the zerocrossing of xk. Equation (31) is reduced to(32)x˙k=fkx,p,u=h2kxi,p+uk≥0.Therefore,xk will increase again as a positive number. As a consequence, negative values have no relevance for the nonlinear model and can simply be corrected in the linearized model. The predictions of the linearized model were thus in the end adjusted by clipping the state variable values in the range 0,∞.


## 2.4. Monte Carlo Analysis of Linearization Accuracy Range

To estimate the quantitative behavior of the linearized model, the following workflow was used:(1)
Generate a random parameter vectorp and set rnoise, representing the range of a uniform distribution centered around 1, used to perturb the input vectors(2)
For each input sampleub,i,(a)
Forward simulate the nonlinear system until the steady-statexss0 is reached(b)
Linearize the nonlinear model aroundxss0(c)
Computeui by perturbing input vector ubf,i, with a multiplicative noise of range ±rnoise%(d)
Forward simulate system (11) until the steady state, using the initial conditions xagreement0=0xss0(3)
Compute metrics by averaging over all input samples, for the currentp and rnoise.


## 3. Results

For the WNT model described Introduction, the proposed workflow was applied four times, for different parameter vectorsp. Three ranges rnoise of multiplicative noise were used: 10%, 25%, and 50%. Vectors p were drawn from a log-uniform distribution over 10−1,101. The same set of inputs ubf,i was used for all workflow runs, and it consisted of 2000 synthetically generated input samples, drawn from a fitted log-normal distribution over a smaller set of real input samples. The perturbed input was computed as(33)ui=1+rnoise−2rnoisedi∘ubf,i,where “∘” is the Hadamard product and di is the perturbation direction for input ui. Vectors di were drawn from a uniform distribution over 0,1. Each workflow run had its own set of directions di, used for all ranges rnoise.First of all, we refer to runtime considerations: a runtime larger than 0.150 seconds is required for a nonlinear simulation to steady state of the herein considered WNT model (based on Matlab’s ODE15 s function, with user-supplied Jacobian), whereas for the linearized system, only 0.003 seconds are required. This time difference may become significant, especially when the “relative” dataset chunk size is considerable.Tables1 and 2 display the steady-state results (averaged over all input samples) for the four runs of the proposed workflow. The predictions output by the linearized model were corrected as described above, to avoid negative-state variable values. For Table 1, the results were averaged over all model states, while for Table 2, the results were averaged over the model outputs (i.e., the readout model states). MRE and MSE represent mean relative and mean squared errors, respectively:(34)MSE=xss linear−xss nonlinear2¯,(35)MRE=xss linear−xss nonlinearxss nonlinear¯,where xss linear and xss nonlinear are the linear and nonlinear steady-state vectors (after having applied perturbation) and stde/stdxdiff is the ratio between the standard deviation of the linearization error vector and the standard deviation of the vector xdiff=xss nonlinear−xss0, where xss0 is the linearization point.Table 1
Metrics on all states, computed using adjusted linear predictions.
Workflow number
1
2
3
4
r
noise
=
±
10
%
std
e
/
std
x
diff (%)
4.61
4.89
5.08
5.53
MRE (%)
0.41
0.27
0.59
0.31
MSE
0.162
0.024
0.195
0.039
r
noise
=
±
25
%
std
e
/
std
x
diff (%)
10.41
12.09
11.60
12.90
MRE (%)
2.33
1.65
3.38
1.86
MSE
4.763
0.895
6.217
1.044
r
noise
=
±
50
%
std
e
/
std
x
diff (%)
18.67
22.68
20.43
22.92
MRE (%)
9.10
6.13
12.29
6.48
MSE
62.042
13.702
85.323
17.656Table 2
Metrics on model outputs, computed using adjusted linear predictions.
Workflow number
1
2
3
4
r
noise
=
±
10
%
std
e
/
std
x
diff (%)
9.65
7.05
6.95
7.33
MRE (%)
0.87
0.40
0.72
0.39
MSE
6.80E−04
1.08E−04
4.50E−05
3.30E−05
r
noise
=
±
25
%
std
e
/
std
x
diff (%)
22.51
17.53
16.76
17.34
MRE (%)
5.26
2.62
4.43
2.61
MSE
2.70E−02
4.52E−03
1.32E−03
1.28E−03
r
noise
=
±
50
%
std
e
/
std
x
diff (%)
43.12
31.61
30.73
31.81
MRE (%)
17.96
10.12
17.80
10.76
MSE
3.89E−01
6.06E−02
2.76E−02
1.95E−02As an example, we display in Figure2 the evolution in time of the error metrics MRE and MSE for a single case with rnoise=50%. The errors increase in time, as the state vector moves farther away from the linearization point xss0, and consequently, the actual Jacobian differs more substantially from the one computed at xss0.Figure 2
Evolution in time of the error metrics for one perturbation example withrnoise=50%.In Figure3, we display the evolution in time for two model states for which the linearization introduces errors (the same case as in Figure 2 was considered). At first, the two models display almost perfect agreement (both qualitatively and quantitatively). However, as the state vector starts departing farther away from the linearization point (at t=0), the error between the two models increases.Figure 3
Evolution in time of two state variables, for both systems, withrnoise=50%.Figure4 displays the dependence between linearization error and rnoise for one state variable of the model. The evolution in time of both models, linearized and nonlinear, is depicted, for all three values of rnoise. For small perturbations, the linearized model follows closely the nonlinear version. As the perturbations become larger in magnitude, the difference between the time courses of the two models increases. This behavior can be observed also in Figure 5, where the MRE variation (for the output variables of the model, in steady state) was plotted for all workflow runs. An approx. polynomial dependence can be observed, which is in agreement with equation (28).Figure 4
Evolution in time of one state variable for the linear and the nonlinear systems, for all three perturbation levels. The differences between the two models are indicated with arrows. The linear model is depicted with darker colors. Forrnoise=10%, the two models overlap perfectly.Figure 5
Dependence between perturbation levelrnoise and the mean relative linearization error of the model outputs, at steady state and for the four workflow runs.Without the adjustment introduced for the predictions of the linearized model, for some input perturbations, certain state variables in the nonlinear model decay to zero, while the state variables in the linearized model converge to negative values. An example is displayed in Figure6. The two models display a good agreement up to the zerocrossing of the state variable in the linearized model. Without the adjustment, this example would have yielded a high MRE value since the linearization error would have been divided by a small denominator, according to equation (35).Figure 6
An example which yields negative values for one of the linearized model state variable. The nonlinear state variable decays to zero, while the state variable in the linearized model converges to a negative value.


## 4. Conclusions

Linearization is a useful tool for the analysis of large nonlinear dynamical systems. It produces a simpler version of the original model, on which established linear analysis methods can be applied to gain insights into the local behavior of the original model.As an example to test the approach, a model of the canonical WNT signaling pathway was considered. The model was forward simulated until steady state was reached for an input setub,i. Linearization was performed at steady state for each input sample, and the accuracy of the linearized model was determined for various ranges of input perturbation magnitudes.The concepts related to linearization validity were presented. A linearization approach offers correct insights only when tested against valid input perturbations (i.e., ones that do not drive the nonlinear system into an unstable region).In terms of accuracy, errors tend to increase in a polynomial fashion w.r.t. the magnitude of tested perturbations. For small enough perturbations, the linearized system closely follows the original one.By choosing an acceptable error level, a relative input perturbation threshold can be derived, which can be used to discriminate between various test input vectors, i.e., predict which input vectors will not lead to large linearization errors. For the considered WNT model, linearization errors were determined to be low for input vectorsui within ±10% of the base input vector ubf,i.



---

*Source: 8461820-2019-06-10.xml*