Stability Constraints of Markov State Kinetic Models Based on Routh-Hurwitz Criterion

In computational neuroscience, receptors, channels and more generally signaling pathways are often modeled with Markov state models to represent biochemical reactions, which are then implemented with bilinear equations. One of the goals of these models, once calibrated with experimental results is to predict the dynamics of the biological system they represent in response to molecular perturbations and therefore facilitate and enhance the success rate of drug discovery and development. To model receptors under both pathological and physiological conditions, modelers usually modify the ligand association and dissociation parameters in the kinetic model during the optimization phase of model development. However, some parameter values may lead to unstable models, (cid:80)(cid:68)(cid:78)(cid:76)(cid:81)(cid:74)(cid:3)(cid:70)(cid:68)(cid:79)(cid:76)(cid:69)(cid:85)(cid:68)(cid:87)(cid:76)(cid:82)(cid:81)(cid:3)(cid:89)(cid:72)(cid:85)(cid:92)(cid:3)(cid:71)(cid:76)(cid:73)(cid:191)(cid:70)(cid:88)(cid:79)(cid:87)(cid:15)(cid:3)(cid:87)(cid:76)(cid:80)(cid:72)(cid:16)(cid:70)(cid:82)(cid:81)(cid:86)(cid:88)(cid:80)(cid:76)(cid:81)(cid:74)(cid:3)(cid:68)(cid:81)(cid:71)(cid:3)(cid:76)(cid:81)(cid:72)(cid:73)(cid:191)(cid:70)(cid:76)(cid:72)(cid:81)(cid:87)(cid:3)(cid:69)(cid:72)(cid:73)(cid:82)(cid:85)(cid:72)(cid:3)(cid:83)(cid:72)(cid:85)(cid:73)(cid:82)(cid:85)(cid:80)(cid:76)(cid:81)(cid:74)(cid:3)(cid:83)(cid:85)(cid:72)(cid:71)(cid:76)(cid:70)(cid:87)(cid:76)(cid:89)(cid:72)(cid:3)(cid:76)(cid:81)(cid:3)(cid:86)(cid:76)(cid:79)(cid:76)(cid:70)(cid:82)(cid:3)(cid:86)(cid:87)(cid:88)(cid:71)(cid:76)(cid:72)(cid:86)(cid:17)(cid:3)(cid:44)(cid:81)(cid:3)(cid:82)(cid:85)(cid:71)(cid:72)(cid:85)(cid:3)


Introduction
One of the main goals of systems biology consists in providing an integrative description of living systems and organisms by using simulation of complex and interconnected mathematical models of metabolic networks and signaling pathways [1,2]. System Biology Markup Language (SBML) [3] constitutes a standard for modeling systems biology and uses an XML-based language. SBML inherently allows developing synaptic receptor models at the molecular level using kinetic parameters as well as integrating a variety of mechanisms, which modulate the physical and functional properties of receptors (i.e., rise time, decay time, amplitude, etc.). Numerical tools are used to simulate these synaptic receptor models.
Rhenovia Pharma has developed a simulation platform (RHENOMS™) for hippocampal glutamatergic synapses and its integration into complex neuronal networks by combining kinetic receptor models and models of other reactions taking place within synapses [4].
is platform has been used to perform in silico experiments [5] to better understand the dynamics of synaptic transmission under physiological and pathological conditions. In order to calibrate these models, it is necessary to optimize the model parameters such that simulated results reproduce selected experimental results. In this optimization phase, the modeler iteratively modi es model parameters until simulated results match experimental ones ( Figure 1). In addition, in order to perform in silico simulations under normal and pathological conditions to test the e ects of drugs for neurodegenerative or other central nervous system (CNS) diseases, both 'physiological' and 'pathological' synaptic receptor models have to be developed. To characterize the behavior of those 'control' and 'diseaselike' models, users go through an initial parameter optimization phase, which is very important to ensure that the model displays the desired behavior under both control and pathological conditions. However, the chosen parameter values may lead to an unstable model. To address this issue, stability conditions could need to be investigated from the mathematical structure of the kinetic model.
From a modeling point of view, biological kinetic models [6][7][8], particularly synaptic receptor models [9][10][11], represent states and transition probabilities between the di erent states (i.e., open or desensitized probability states). e resulting model, called Markov kinetic model, is o en a ne in control, which is a particular form of nonlinearity. ese kinds of models are known as bilinear systems and Studying the kinetic model to nd the proof of stability is not enough in this presented model tuning process. During the model tuning process, users modify the model associations and dissociations parameters (with an optimization algorithm). e problematic encountered is that some parameter values and simulation conditions may lead to model instability resulting in (i) lack of convergence of the simulation runs, (ii) erroneous solutions of the sets of ODE used to solve the model equations, and (iii) prohibitively long simulation durations.
is instability may jeopardize model simulation and consequently hinder the success of the parameter optimization phase. ese issues led us to ask whether it could be possible to identify constraints on parameter values, which could guarantee model stability. Finding these stability constraints is the primary objective of this study in order to guarantee the optimization step success.
In this manuscript, the stability concept and its application for linear and bilinear models are presented in Section 2.
en, the proposed approach based on the Routh-Hurwitz stability criterion is presented. e RH criterion e cacy is demonstrated with two examples, a textbook model and the GABA A receptor model. Finally, discussion about the results obtained with the RH criterion as well as their limitations and elaborates on the perspectives of this work.

Model Stability
To represent the notions of stability and equilibrium point one must rst de ne the system of interest. Let's consider a process represented by the following nonlinear equation: If an equilibrium point e x is attractive and stable, then it is said to be asymptotically stable.
Otherwise, the equilibrium point e x is considered to be unstable [24]. x t

Stability Studies for Linear System
x [24]. In other words, a system is said to be stable if it returns to its steady state a er a disturbance by an are commonly modeled by an Ordinary Di erential Equations (ODE) [12][13][14] as follows: where ( ) A is a real constant matrix. Notice that the analytical solution of this system cannot be o en obtained due to the complexity and nonlinearities of the system. However under some conditions, stability conditions are proposed, as mentioned by Chen et al. [15]. Most studies consider time-invariant continuous bilinear systems with linear feedback [16][17][18]. Almost all of these suggest studying stability by nding a su cient condition for the existence of a feedback control, such that the resulting closed-loop system is asymptotically stable.
A large number of publications study the stability of bilinear systems with the Lyapunov functions [19,20].
is energy based method developed by Lyapunov in 1892 [21], consists in nding su cient conditions for stability of nonlinear systems. e su cient conditions require the existence of Lyapunov functions ( ) u t for the system or the existence of some negative (or positive) de nite matrices, such that the feedback input ( ) u t depends on those matrices as shown in ref. [22,23]. is method avoids the explicit solution of the ODE and can then be applied to a large class of nonlinear systems. is method could represent an interesting way to solve the problem investigated in this paper. However, for each modi cation of a given parameter, another function ( ) V x needs to be determined during the optimization phase, increasing the computation-time. In order to reduce the mathematical complexity of the investigated model, an approximation (i.e., a linearization) around an operating point can be proposed considering the system input ( ) u t as a piecewise constant signal. In this way the nonlinear model is transformed into a linear system and some well-known stability conditions can be easily employed based on the computation of the eigenvalues of the dynamic matrix.
Among them, Routh-Hurwitz (RH) stability criterion is a mathematical test, which provides a necessary and su cient condition to characterize the stability of a linear system without computing explicitly the solution of the eigenvalues (roots), thus reducing computing time. is test is well adapted to biological bilinear systems, which are considered here because they may become extremely complex, consequently making it di cult (computationally speaking) and sometimes even impossible to compute the eigenvalues (roots). In addition, it is also important to note that even if the proof of stability modeling, generates an initial, not-validated model. To obtain a physiologic or pathologic model, users modify the model parameters during the optimization step. The optimization step is run until desired physiological or pathological model is reached. The third step, or the study step, represents the application stage of the model, when new predictions can be generated by the model. column shows the existence of a new positive real part for the roots of the system studied. e appearance of a zero in the pivot column indicates the existence of a pure imaginary root [25,26]. In that case, the studied system corresponds to an oscillating system also called pseudo-stable. e Routh-Hurwitz criterion presented in this section represents a method of determining the location of the roots of a polynomial with constant real coe cients with respect to the le and right halves of the complex plane (also called s-plane), without actually solving for the roots.

Stability Constraints for a Kinetic Model
e contribution work, presented in this section, is inspired from RH criterion to nd stability constraints of a linearized kinetic model. To extract constraints guarantying the kinetic model stability, the following algorithm is proposed: Linearize the system (1) around an operating point considering the input as a piecewise constant 2 Compute the characteristic polynomial (6) of the linearized system 3 Construct the RH table with equation (6) and rules (7)  4 e rst column (pivot column) of the table contains the constraints needed to guarantee the model stability e kinetic model will be stable if all constraints provided by the RH criterion have the same sign (positive or negative). Remark that the constraints will provide only a local stability (around the chosen operating point) due to the linearization.
e following examples illustrate the proposed algorithm for designing stability constraints.

Textbook example
To illustrate the RH criterion for nding stability constraints on kinetic model parameters, a very basic kinetic model as shown in Figure  3 is used. is kinetic scheme represents a textbook case composed of three states 1  x k x t k x t u t x k x t u t k x t k k x t x k x t k x t (8) with k 1 =0.7, k 2 =0.2 and k 3 =0.4 outside element. Notice however that establishing stability conditions for nonlinear systems is an arduous task. O en, a linear approximation of the nonlinear system is employed in order to verify local stability of the nonlinear system around an operating point. In this paper, the kinetic model input is considered as a piecewise constant signal, which amounts to linearizing the system around an operating point.
A linear system (where A is the dynamic matrix), for which the analytical solution is is negative. More explicitly, for a linear system, the real part of the eigenvalues ofA ( ( 0 , where represents the eigenvalues), needs to be negative to ensure stability.
It can be recalled here that nding the eigenvalues of the matrix A, with det I A , is equivalent to nding the roots of the polynomial characteristics ( ) P de ned by: Where λ represents the A matrix eigenvalues and a i the polynomial characteristic constant coe cients. e characteristic polynomial degree will be the same as the kinetic model size, i.e., the number of the states. Based on the complexity of the linear system, computation of the roots of the polynomial characteristics can be di cult or impossible. In the following section, an interesting stability test is presented to avoid the explicit computation of the roots.

Routh-Hutwitz Criterion
In the eld of automatic control, the Routh-Hurwitz stability criterion is a mathematical test, which provides a necessary and su cient condition to characterize the stability of a linear system without computing explicitly the roots of (6). e Routh-Hurwitz test is a recursive algorithm to determine if all the roots of (6) have negative real parts [25]. e Hurwitz Matrix (or Table) helps to conclude on the stability of the linear system: the system is said stable if and only if the determinants of its principal submatrices are all positive [25]. e crux of the Routh-Hurwitz criterion [25] lies in the roots λ of the characteristic polynomial of a linear system. ese roots with negative real parts represent the stable solutions ( ) t e of the system (bounded solutions). e exact calculation of the roots is then replaced by the Routh-Hurwitz criterion whether the roots are positive or negative. us, the test provides a means to determine whether the steady state point of a linear system is stable avoiding to directly solving the di erential equations that characterize it. e rst two lines of the RH table are lled in column, with the coe cients of the characteristic polynomial (6). e rst row contains the coe cients of the terms in Note that an empty cell is considered equal to zero, if necessary. e calculation of the various cells of the table continues until the rst column is completed. e rst column of this table is called the pivot column.
e Routh-Hurwitz criterion helps concluding whether a system is asymptotically stable when all elements of the pivot column have the same sign. Each change of sign in the elements of the pivot  (9) is: e Routh-Hurwitz table of the simpli ed model is built with the n a parameters of equation (10) and is presented in Table 1 with 3 1, 3 : 0 cst k u k k k k (13) e purely imaginary root of this textbook example, represented by c 0 does not provide an asymptotic stability for this model (see Section 2). Only a simple stability can be guaranteed for the three-state system. For example by choosing the parameters k 1 =0.7, k 2 =0.2, k 3 =0.4 and k 4 =0.4 0 =1.0 [28], the constraints necessary to have a stable model are met and a simulation around this operating point is shown in Figure 4. An input 0 =1.0 is applied at 50 msec during 20 msec, and is shown in Figure 4A (black). During this interval (50-70 msec), the states respectively presented in Figure 4B, 4C and 4D (red, green and blue) change the value of concentration without generating divergence. Once the application of the input u 0 is interrupted, each state returns to its respective steady state.
It should be noted that some models might receive negative input signals, such as voltage-dependent channels. For this reason, the three-state model was also tested while retaining 0 parameters and applying a negative input value equal to -0.5. With this value of 0 it is easy to see that the stability constraints speci ed by equations (12) and (13) are not met; therefore the model is unstable for a parameter 0 =-0.5. e instability of the model as provided by the Routh-Hurwitz criterion is veri ed by the simulation presented in Figure 5. Indeed, x 1 , x 2 and x 3 represent the model states and an input to the system. The k i parameters represent transition probabilities between the different states.     Figure 5C and 5D of the system states is not possible [28].
When the system is no longer excited by the negative value of 0 , each state returns to its respective steady state. e stability constraints (equations (11) to (13)) provide a signi cant advantage for parameter optimization. Indeed, the constraints provided by the RH criterion allow determining the range of validity of a varying parameter. It is therefore possible to determine an area of stability in the studied neighborhood. In the case of the linearized three-state model (equation (9)), the validity domain of the parameter 0 ensuring stability of the system which depends on other parameters can be obtained by replacing the 0 2 3 0 0.06 u k k u parameters values in equations (11) to (13). Easily, it may conclude that to ensure the stability of this model, the parameter 0 must satisfy the following inequalities:

GABA A receptor model
A er this textbook case example, the proposed algorithm is applied to a GABA A receptor model. e GABA A receptor is an ionic receptor and is activated by application of gamma-aminobutyric acid, aka GABA. is ionotropic receptor is of great importance since it binds the main inhibitory neurotransmitter in the CNS and constitutes an important target in the treatment of many CNS diseases [29,30]. Figure 7 shows the kinetic scheme of the GABA A receptor model developed [30] and used for this study. is model consists of eight states represented by rectangles on the kinetic scheme. A parameter optimization phase is considered to modify the GABA A receptor model deactivation to perform the stability analysis. In other word, optimization of association and dissociation parameters outlined in red in Figure 7 will be conducted to determine stability constraints according to these parameters.
Considering the GABA input as a piecewise constant GABA 0 , it is possible to linearize the system and therefore obtain the following matrix: By arranging the parameters of (15) in Table 2 and by applying the Routh-Hurwitz criterion on this table, the constraints that ensure the stability of the model in the rst column of this table (pivots column) are obtained. e stability domain of the GABA A receptor model is de ned by the pivot column of the Routh-Hurwiz table (Table 2). e parameters in this column represent the constraints to be respected to guarantee the stability of the GABA A receptor model during the parameter optimization phase. e parameter a 8 is equal to 1, implying that all constraints must be positive to ensure stability of the equilibrium point of the model. However, s parameter is equal to zero; which implies that an asymptotic stability cannot be provided. Only stability around the operating point may be guaranteed for this model.
To achieve a deactivation optimization on GABA A receptor model ka 3 , kd 3 , ka 5 , ka 8 and kd 8 parameter values need to be modi ed. ey represent the optimization settings for our stability study. Using the parameter values given by ref. [31], the stability constraints provided by the Routh-Hurwitz criterion are met and the response of the simulation remains stable for 1mM of GABA during 250 msec as   Figure 8. e optimization parameters are changed randomly between 0 and +30. is range of values corresponds to the domain of validity of the association and dissociation rate constants proposed by ref. [31] for our GABA A model receptor. is random change in optimization settings places the system in situations where the stability constraints are no longer guaranteed and will generate unstable simulations conditions. Figure 9 illustrates responses of the model in an instable condition, when 1mM of GABA is applied during 20 msec. Instability on the desensitization parameters of the GABA A model receptor produces instability on desensitization state 'DA2F'. As shown in Figure 9D, instability on 'OA2' open state appears too, suggesting that the instability propagates in the kinetic model. 'DA2F' state is biologically unstable, as a negative concentration is impossible; 'OA2' state is unstable too: it corresponds to the channel opening probability for the receptor and a probability superior to 1 is therefore impossible (1 corresponds to 100% opening for the receptor).
To perform an analysis of the GABA A receptor model, the RH criterion enables the identi cation of stability constraints during parameter optimization phase. Indeed, once the model is linearized around the operating point of interest to the user, it is possible to identify constraints that ensure model stability.
e constraints may then be integrated within the optimization process making the operation completely transparent to users regardless of their area of expertise.
Besides determining model stability in the parameter optimization phase this criterion may also be used to determine the stability domain of model parameters found in the literature. Indeed, keeping the GABA A model receptors with the parameters used in ref. [31], it is possible to determine whether the supplied parameters allow the model to operate in stable conditions. Figure 10 shows the stability constraints of the GABA A receptor model with parameters used in ref. [31]. To ensure the stability of the model all constraints must have the same sign. Maintaining the constraint a 8 to a constant equal to 1 therefore implies that all other constraints must remain positive. us, the   The parameters used were obtained after random change between 0 and 30 of the parameters considered during desensitization optimization. 1mM GABA is applied as input to the system for 20 msec. The associations and dissociations parameters are provided from [31]. All constraints are positive or zero between the two red dotted lines, indicating that with [GABA] between 0 and 6.54mM, the GABA A receptor model is stable.  Figure 10 shows that all constraints are positive or equal to zero between the two dotted lines. Indeed, h 0 constraint is a constant equal to zero, implying that this GABA A receptor model may not be asymptotically stable, but globally stable only if all other constraints are positive. e dotted line on the le corresponds to a GABA concentration equal to -0.0015 mM while the one on the right corresponds to a GABA concentration equal to 6.54 mM. Indeed, for a numerical value of [GABA] less than -0.0015 mM, the constraints d 4 , e 3 and f 2 are negative; whereas for a numerical value of [GABA] greater than 6.54 mM the constraint f 2 is again negative. It is obvious that no negative GABA concentration may be administered. Finally, the stability domain reading s of this GABA A receptor model with the parameters proposed in ref. [31] is straightforward, as described in Figure 10. e stability domain s of our studied bilinear kinetic receptor model can be written as follows: : 0 ; 6.54 s GABA (16) e Routh-Hurwitz criterion therefore leads us to conclude that for the GABA A receptor model with the parameters proposed by Pugh and Raman [31]. e model is globally stable if the concentration of GABA remains between zero and 6.54 mM.

Conclusions
e primary objective of this study was to determine stability constraint to facilitate the development of kinetic models. We proposed to linearize Markov state kinetic models around an operating point considering the input of the bilinear model as piecewise. is linearization allowed application of the Routh-Hurwitz criterion. We rst applied this criterion on a textbook case model and demonstrated that the resulting constraints indeed guaranteed stability of the model, thus increasing the chances of success of the optimization phase. Similarly, without the constraints equations guaranteeing the stability, a random parameters change in optimization process may result in instability in the model. In addition to identifying constraints that ensure the stability of synaptic receptor models, the Routh-Hurwitz criterion also allows identi cation of the stability range, i.e., the determination of the range of input values allowed for a prede ned parameter set. We illustrated our approach with a brief study of a GABA A receptor model and showed that the parameters provided by Pugh and Raman [31] de ne a model that is numerically stable if GABA concentration remains under 6.54 mM.
It is important to remember that bilinear kinetic receptors are modeled by ODE. e resolution of an ODE is discretized; this implies that bilinear kinetic receptor models are discrete time systems. As the studied models are continuous but resolved by a discrete time method then the stability constraints needs to be obtained with a discrete time criterion. e Jury criterion [32] is the equivalent of RH criterion for linear discrete time systems. Remark that bilinear kinetic receptor models are mostly hybrid [33]. Hybrid systems combine several operation modes. Voltage-dependent calcium channels (VDCC) are an example of a well-known hybrid system in neuroscience [34]. As presented by Ja e et al. [34], the equations describing the VDCC model generate a speci c behavior for a membrane potential that is negative or equal to zero and a completely di erent behavior when the potential becomes strictly positive. Future work will consist in adapting the RH and Jury criteria methods to such bilinear hybrid systems. Notice that computational neuroscience o en involves models of ever increasing computational complexity [35], our team has already invested signi cant e orts in reduction of model complexity. Another orientation work may be considered, in parallel to the studied stability constraints here, in investigating a reduction of model order [36,37], while guaranteeing its stability and maintaining desired nonlinear dynamic response.