#### 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.

#### 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.

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 as , 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:where is the state vector, is the parameter vector, is the input vector, and is the output vector. Function describes the behavior of the system, while function defines the output variables of the model. If is Lipschitz continuous, then equation (1) has a unique solution [9].

For the considered WNT model, and are Lipschitz continuous functions and parameter vector is presumed constant (therefore, its notation can be omitted). The model has the following size: , , , and .

If input vector has a step-like evolution, then at a fixed point (considered as steady state), equations (1) and (2) are rewritten as follows:where is the steady-state output vector and is the final value of the base input vector .

If a linearization is performed around , using Taylor series expansion, and all terms of second or higher order are being discarded, thenwhere is the state Jacobian evaluated at and is the input Jacobian evaluated at . Vectors and 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]:

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 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 on or )(ii)The state Jacobian does not depend explicitly on

Using equations (3) and (6)-(7),where is the state vector of the linearized system.

Based on equation (3), at the point of linearization, is zero. Therefore,

If is defined as , i.e., the vector representing the deviations of the state from the point of linearization for the linear model, then

The agreement between the linearized model and the original model can be formalized as follows:where is the state evolution as dictated by the original nonlinear system and 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 definingan analysis related to the errors introduced by the linearization can be conducted for different input vectors . An example workflow would be as follows:(a)Forward simulate the nonlinear model until a steady-state point is reached for input (b)Compute matrices and (c)Forward simulate model (11) for a new input vector , with the initial conditions (d)Compute measures of interest on the agreement error vector

##### 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 is “hyperbolic” if all eigenvalues of the Jacobian evaluated at (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:

Herein, all vector fields which have the following properties are regarded as “step-like” inputs:

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 conditions .

However, these considerations are not universally true for nonlinear systems. The existence of the steady state is guaranteed if the following lemma applies. Suppose is a nonempty bounded subset of . By definition [13], the limit of set , denoted , is the set of all points for which a sequence of pair exists, with and such that

Lemma. *[13]. Let be a nonempty bounded subset of and suppose there is a number such that for all and all . Then, is a nonempty compact invariant set. Moreover, the distance of from approaches 0 as , uniformly in .*

Less formally stated, if the nonlinear system has an overall stable behavior, i.e., with initial conditions in a neighborhood , its state does not diverge towards infinity, then at a large enough time , its state trajectory will enter and remain in a closed set . If neighborhood is centered around an attractor fixed point , 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 consist of only one point, i.e., the system has no oscillatory behavior at steady state.

To analyze the influence of the initial condition 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 , with the property that , where and 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 with respect to direction is defined as the limit (if it exists) [15]:where . The global Lyapunov exponent describes the decay in time of the distance between two trajectories which start at two separate initial points and :

If , any state trajectory which starts at (for any sufficiently small ) will converge to the “main” trajectory which starts at . Since for a small we can state thatwhere is the Jacobian of evaluated at and can be expressed as

Thus, the global Lyapunov exponent of at point with respect to direction is the average of the local Lyapunov exponents of at points of the trajectory of w.r.t direction at each point [15]. The local Lyapunov exponents are characterized by the eigenvalues of . Therefore, if the system reaches a new steady state (e.g., an attractor fixed point ) when an input perturbation is applied (starting from an initial stable equilibrium point which can be employed for the linearization), all state trajectories starting from a sufficiently small neighborhood of will converge to .

Let be a compact and convex set of points for which the state Jacobian evaluated at has negative real-part eigenvalues. More formally, let be the eigenvalues of , and . In this case, the eigenvalues of the discrete-time model lie inside the unit circle. is restricted so that . Then, at each point , the local Lyapunov exponent w.r.t. direction is defined asand is negative; therefore,where and is a small enough number.

It can be shown that, for a constant input , for every initial point , the nonlinear model reaches a unique fixed point at steady state. Banach’s fixed point theorem states that if is a strict contraction, then has a unique fixed point in [9]. A map is a strict contraction on , if :where . To prove equation (22), consider the line segment connecting points and , and the sequence of distinct points between and . Then,

The points of each pair (, ) are close enough:where is the local Lyapunov exponent of at . Applying the procedure of equation (23) iteratively and employing equation (24) yields

Let .

Then,

Because are collinear points, their inner distances sum to . Therefore,and is a strict contraction for all . This also implies that the limit of consists of only one point. For comparison, a stable linear system will always reach a unique steady-state fixed point (i.e., at ) 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 is a particular case of the global Lyapunov exponent reasoning described above. may also be regarded as a smaller subset of the attraction basin of the attractor fixed point .

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 point , 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 , before applying the exogenous perturbation). Although the trajectories of the models will initially behave qualitatively the same when starting from , i.e., the trajectories will start varying in the same direction (since the linearization is accurate for small deviations around ), the two models may behave qualitatively differently around . This holds true especially if the Jacobian varies considerably between and .

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, let be an attractor fixed point around which a linearization was performed, and the input perturbation be a Heaviside step function. As a first example, consider the Taylor expansion of around :

At a new steady state of the nonlinear model , 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 w.r.t. at both fixed points and yields

If the Jacobian varies significantly between the two fixed points and , the linearized and the original system will exhibit different input sensitivities. The rate of change of the Jacobian w.r.t. 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 so that, for every , if , then . At the end of each time increment, the linearization (i.e., the Jacobian) is recomputed. Applying this for steps until the new steady state is reached yields

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 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 function has for all state variables the following form:where is the *k*^{th} state variable, is the set of all state variables excluding (i.e., ), and and are nonlinear Lipschitz continuous functions with the property that and , , and is constant. When both the exogenous input and the parameter vector are chosen to be positive (which is the case for the considered datasets) and all state variables are nonnegative, the subsystem defined by 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 (at ), will always remain nonnegative. To prove this, let all other state variables be nonnegative and consider the zerocrossing of . Equation (31) is reduced to

Therefore, 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 .

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

#### 3. Results

For the WNT model described Introduction, the proposed workflow was applied four times, for different parameter vectors . Three ranges of multiplicative noise were used: 10%, 25%, and 50%. Vectors were drawn from a log-uniform distribution over . The same set of inputs 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 aswhere “” is the Hadamard product and is the perturbation direction for input . Vectors were drawn from a uniform distribution over . Each workflow run had its own set of directions , used for all ranges .

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.

Tables 1 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:where and are the linear and nonlinear steady-state vectors (after having applied perturbation) and is the ratio between the standard deviation of the linearization error vector and the standard deviation of the vector , where 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 . The errors increase in time, as the state vector moves farther away from the linearization point , and consequently, the actual Jacobian differs more substantially from the one computed at .

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 ), the error between the two models increases.

Figure 4 displays the dependence between linearization error and for one state variable of the model. The evolution in time of both models, linearized and nonlinear, is depicted, for all three values of . 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).

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. 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).

#### 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 set . 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 vectors within of the base input vector .

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CCCDI—UEFISCDI, project number ERANET-FLAG-ITFoC (2), within PNCDI III.