Effect of Linearization in a WNT Signaling Model

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.


Introduction
e WNT signaling pathway is a cascade of biochemical reactions transducing signals from the extracellular space into the cells. is pathway is involved in different biological processes such as embryonic development, tissue homeostasis, and tumorigenesis [1][2][3]. e 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. e last two components phosphorylate β-catenin on several serine and threonine residues targeting it for ubiquitination and subsequent proteosomal degradation. is continual elimination of the β-catenin prevents its translocation to the nucleus [1][2][3]. e pathway is activated when WNT ligands bind FZD and LRP5/6. is 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.
e 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]. e 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. e 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. ese 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. is 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).
is 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][7][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. ese 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 as u b , and the other input vectors are regarded as input perturbations.

Nonlinear Systems and Linearization.
Generically, a dynamic nonlinear system is described using the following equations: 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). e model has the following size: x ∈ R 421 , u ∈ R 60 , p ∈ R 1154 , and y ∈ R 40 .
If input vector u has a step-like evolution, then at a fixed point x ss (considered as steady state), equations (1) and (2) are rewritten as follows: where y ss is the steady-state output vector and u bf is the final value of the base input vector u b . If a linearization is performed around x ss , using Taylor series expansion, and all terms of second or higher order are being discarded, then where (zf/zx)x ss , u bf is the state Jacobian evaluated at (x ss , u bf ) and (zf/zu)x ss , u bf is the input Jacobian evaluated at (x ss , u bf ). Vectors (x − x ss ) and (u − u bf ) 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]: 2 Computational and Mathematical Methods in Medicine 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 of (u − u bf ) since they do not exist. ese conditions are true for the WNT model considered in this paper: (i) e input Jacobian is constant (i.e., does not depend on x or u) (ii) e state Jacobian does not depend explicitly on u Using equations (3) and (6)- (7), where x linear is the state vector of the linearized system. Based on equation (3), at the point of linearization, dx ss /dt is zero. erefore, d dt If x dev is defined as (x linear − x ss ), i.e., the vector representing the deviations of the state from the point of linearization for the linear model, then e agreement between the linearized model and the original model can be formalized as follows: where dx nonlin /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 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:

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 x f is "hyperbolic" if all eigenvalues of the Jacobian evaluated at x f (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. Computational and Mathematical Methods in Medicine 3 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: Herein, all vector fields u which have the following properties are regarded as "step-like" inputs: is 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 conditions x 0 .
However, these considerations are not universally true for nonlinear systems. e existence of the steady state is guaranteed if the following lemma applies. Suppose B is a nonempty bounded subset of R n . By definition [13], the ω limit of set B, denoted ω(B), is the set of all points x for which a sequence of pair x k , t k exists, with x k ∈ B and lim Lemma. [13]. Let B be a nonempty bounded subset of R n and suppose there is a number M such that Less formally stated, if the nonlinear system has an overall stable behavior, i.e., with initial conditions x 0 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 x ss , 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 condition x 0 on the state trajectory of the nonlinear model, the Lyapunov exponents theory can be employed. We consider Φ as the discretetime counterpart of the continuous-time model f, with the property that x k+1 � Φ(x k ), where x k � x(t 0 + kT s ) and T s is the sampling period. e 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]: where Φ k (x 0 ) � x(t k , x 0 ). e 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: 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 where J(Φ, x) is the Jacobian of Φ evaluated at x and λ can be expressed as us, the global Lyapunov exponent of Φ at point x with respect to direction y is the average of the local Lyapunov exponents of Φ at points x 0, x 1 , x 2 , . . . of the trajectory of x w.r.t direction J(Φ k , x)y at each point x k [15]. e local Lyapunov exponents are characterized by the eigenvalues of J(Φ, x k ). erefore, if the system reaches a new steady state (e.g., an attractor fixed point x ssfinal ) when an input perturbation is applied (starting from an initial stable equilibrium point x ss0 which can be employed for the linearization), all state trajectories starting from a sufficiently small neighborhood B of x ss0 will converge to x ssfinal .
Let S ⊆ R 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 (zf/zx)x, p, and R(λ 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. en, at each point x ∈ S, the local Lyapunov exponent w.r.t. direction y is defined as and is negative; therefore, where 0 < c < 1 and Δ > 0 is a small enough number. It can be shown that, for a constant input u, for every initial point x 0 ∈ S, the nonlinear model reaches a unique fixed point x ss ∈ 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: where c < 1. To prove equation (22), consider the line segment connecting points p and q, and the sequence of distinct points p 1 , p 2 , . . . p k , q k , q k−1 , . . . , q 1 between p and q. en, e points of each pair (p k−1 , p k ) are close enough: where c p k < 1 is the local Lyapunov exponent of Φ at p k . Applying the procedure of equation (23) iteratively and employing equation (24) yields Let c max � max c p 1 , . . . , c p k−1 , . . . , c q 1 < 1. en, Because p, p 1 , p 2 , . . . p k , q k , q k−1 , . . . , q 1 , q are collinear points, their inner distances sum to ‖p − q‖. erefore, and Φ is a strict contraction for all p, q ∈ S. is 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., ∀x 0 ∈ R n at t 0 ) 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 of S 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 x ss .
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. e application of an exogenous input perturbation modifies the location in the state space of the attractor fixed point x ss , 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 x ss0 , before applying the exogenous perturbation). Although the trajectories of the models will initially behave qualitatively the same when starting from x ss0 , i.e., the trajectories will start varying in the same direction (since the linearization is accurate for small deviations around x ss0 ), the two models may behave qualitatively differently around x ss . is holds true especially if the Jacobian varies considerably between x ss and x ss0 . e 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, let x lin 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 x lin : At a new steady state of the nonlinear model x ss , f(x ss , 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 of x w.r.t. u at both fixed points x lin and x ss yields If the Jacobian varies significantly between the two fixed points x lin and x ss , the linearized and the original system will exhibit different input sensitivities. e rate of change of the Jacobian w.r.t. δx ≈ x ss − x lin 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 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 If the Jacobian is constant, equation (30) reduces to equation (13). erefore, 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 x lin has a relatively large norm, linearization errors are to be expected.

Properties of the Considered WNT Model
. By analyzing the model equations, we observe that function f has for all state variables the following form:

Computational and Mathematical Methods in Medicine
where x k is the k th state variable, x i is the set of all state variables excluding x k (i.e., i ∈ [1, . . . , k − 1, k + 1, . . . , 421]), and h 1 k and h 2 k are nonlinear Lipschitz continuous functions with the property that h 1 k (x, p) > c k > 0 and h 2 k (x, p) ≥ 0, ∀x, p > 0, and c k is constant. When both the exogenous input u k 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 f k is equivalent to a stable linear time-variant first-order filter.
is, 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 (at t 0 ), x k will always remain nonnegative. To prove this, let all other state variables be nonnegative and consider the zerocrossing of x k . Equation (31) is reduced to erefore, x k 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. e predictions of the linearized model were thus in the end adjusted by clipping the state variable values in the range [0, ∞).

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 vector p and set r noise , representing the range of a uniform distribution centered around 1, used to perturb the input vectors (2) For each input sample u b,i , (a) Forward simulate the nonlinear system until the steady-state x ss0 is reached (b) Linearize the nonlinear model around x ss0 (c) Compute u i by perturbing input vector u bf,i , with a multiplicative noise of range ± r noise % (d) Forward simulate system (11) until the steady state, using the initial conditions x agreement0 � 0 x ss0 (3) Compute metrics by averaging over all input samples, for the current p and r noise .

Results
For the WNT model described Introduction, the proposed workflow was applied four times, for different parameter vectors p. ree ranges r noise of multiplicative noise were used: 10%, 25%, and 50%. Vectors p were drawn from a loguniform distribution over [10 −1 , 10 1 ]. e same set of inputs u bf,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. e perturbed input was computed as where "∘" is the Hadamard product and d i is the perturbation direction for input u i . Vectors d i were drawn from a uniform distribution over [0, 1]. Each workflow run had its own set of directions d i , used for all ranges r noise . 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. is time difference may become significant, especially when the "relative" dataset chunk size is considerable.
Tables 1 and 2 display the steady-state results (averaged over all input samples) for the four runs of the proposed workflow. e 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: MRE � x ss linear − x ss nonlinear x ss nonlinear , where x ss linear and x ss nonlinear are the linear and nonlinear steady-state vectors (after having applied perturbation) and std(e)/std(x diff ) is the ratio between the standard deviation of the linearization error vector and the standard deviation of the vector x diff � x ss nonlinear − x ss0 , where x ss0 is the linearization point.
As an example, we display in Figure 2 the evolution in time of the error metrics MRE and MSE for a single case with r noise � 50%. e errors increase in time, as the state vector moves farther away from the linearization point x ss0 , and consequently, the actual Jacobian differs more substantially from the one computed at x ss0 .
In Figure 3, 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 4 displays the dependence between linearization error and r noise for one state variable of the model. e evolution in time of both models, linearized and nonlinear, is depicted, for all three values of r noise . 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. is 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).
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 Figure 6. e 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).

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. e model was forward simulated until steady state was reached for an input set u b,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. e 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 vectors u i within ±10% of the base input vector u bf,i .

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.