Interval Constraint Satisfaction and Optimization for Biological Homeostasis and Multistationarity

Abstract

Homeostasis occurs in a biological system when some output variable remains approximately constant as one or several input parameters change over some intervals. When the variable is exactly constant, one talks about absolute concentration robustness (ACR). A dual and equally important property is multistationarity, which means that the system has multiple steady states and possible outputs, at constant parameters. We propose a new computational method based on interval techniques to find species in biochemical systems that verify homeostasis, and a similar method for testing multistationarity. We test homeostasis, ACR and multistationarity on a large collection of biochemical models from the Biomodels and DOCSS databases. The codes used in this paper are publicly available at: https://github.com/Glawal/IbexHomeo.

Appendices

Appendix A1: Homeostasis of Low Complexity Models

BIOMD614 is a one species model, with equation:

$$\begin{aligned} \dot{x}=k_1+k_2k_3x-k_1x-k_2k_3x^2 \end{aligned}$$

(7)

At steady state, this leads to:

$$\begin{aligned} k_1+k_2k_3x=(k_1+k_2k_3x)x \end{aligned}$$

(8)

If $k_1\ne 0$, the only solution to (8) is $x=1$, which is the answer given by IbexHomeo. The model describes the irreversible reaction kinetics of the conformational transition of a human hormone, where x is the fraction of molecules having undergone the transition, which is inevitably equal to one at the steady state [14].

BIOMD629 has 2 reactions and 5 species, provides a 2-homeostasis for kinetics parameters with conserved total amounts fixed, provided by the SBML file. This model does not provide ACR, and the homeostasis found can be explained by the conserved total amounts, that lock species to a small interval. But if we change these total amounts and try again an homeostasis test, it should fail. Indeed this model is given by the equations :

$$\begin{aligned} \dot{x_1}=-k_2x_1x_3+k_3x_2, \, \dot{x_2}=k_2x_1x_3-k_3x_2-k_4x_2x_4+k_5x_5, \nonumber \\ \dot{x_3}=-k_2x_1x_3+k_3x_2, \, \dot{x_4}=-k_4x_2x_4+k_5x_5, \, \dot{x_5}=k_4x_2x_4-k_5x_5, \nonumber \\ x_4+x_5=k_6, \, x_2+x_5+x_3=k_7, \, x_2+x_5+x_1=k_8 \end{aligned}$$

(9)

Here $x_3$ (receptor), and $x_4$ (coactivator) have been found homeostatic w.r.t. variations of the kinetics parameters. The total amounts are $k_6=30$, $k_7=\frac{7}{2000}$, $k_8=\frac{1}{2000}$. With these values, we get $x_5 \in \ ]0,\frac{1}{2000}[$, which implies $x_4\in \ ]29.9995,30.0005[$. In the same way we have $x_2+x_5\in \ ]0,\frac{1}{2000}[$, which implies $x_3\in \ ]\frac{6}{2000},\frac{7}{2000}[$. If $k_6,k_7,k_8$ were closer to each other, there would be no reason for homeostasis. This example corresponds to the homeostasis mechanism known in biochemistry as buffering: a buffer is a molecule in much larger amounts than its interactors and whose concentration is practically constant.

Appendix A2: Multistationarity Statistics

Among the 297 models tested for multistationarity using IbexSolve, 63 provide a timeout. For the solved models, 35 do not have steady-state, 153 have a unique steady state, 42 provide multistationarity, 4 have a continuum of steady states. Theses results are given by IbexSolve, but some multistationarity models, such as 003, have in reality an oscillatory behavior (Table 4, 5, 6, 7, 8, 9 and 10).

Table 4. Time statistics for multistationarity

Multistationnary models :

• 10 steady states: 703.

• 7 steady states: 435.

• 4 steady states: 228, 249, 294, 517, 518, 663, 709.

• 3 steady states: 003, 008, 026, 027, 029, 069, 116, 166, 204, 233, 257, 296, 447, 519, 573, 625, 630, 687, 707, 708, 714, 729.

• 2 steady states: 079, 100, 156, 230, 315, 545, 552, 553, 688, 713, 716.

Appendix A3: Homeostasis Results Tables

Table 5. Models with less than 9 species tested for homeostasis w.r.t. the volume of compartments. Minimum and maximum volume of each compartment is set as value given by the SBML file divided by 10 and multiplied by 10, respectively. Kinetic rates and total amount of conservation laws stay fixed. In BIOMD355, all species are 2-homeostatic (and seems independent of the volume). In BIOMD433, the species MK_P is 2-homeostatic.

Table 6. Models with less than 11 species tested for ACR, the second number indicate the number of species that verify 2-homeostasis (ACR included). The minimum and maximum value of each total amount of conservation laws is the value computed from the SBML file, divided by 10 and multiplied by 10, respectively. Kinetic rates and volume compartments stay fixed. FeinbergShinar models serve as tests to confirm that we detect ACR. In BIOMD041, ATP and ATPi are detected to be 2-homeostatic despite the timeout. In BIOMD622, R1B, R1Bubd, and Z are also detected to be 2-homeostatic despite the timeout. In BIOMD413, auxin verify ACR. In BIOMD738, FeDuo, FeRBC, FeSpleen, FeLiver, Hepcidin, FeRest, FeBM are ACR.

Table 7. Models with less than 9 species tested for homeostasis w.r.t. kinetics rates. The minimum and maximum value of each kinetic rate is given by the SBML file divided by 100 and multiplied by 100, respectively. In BIOMD614, the unique species is 2-homeostatic (moreover independent). In BIOMD629, receptor and coactivator are 2-homeostatic.

Table 8. Models with more than 10 species tested for homeostasis w.r.t. the volume of compartments. Minimum and maximum volume of each compartment is set as value given by the SBML file divided by 10 and multiplied by 10, respectively. Kinetic rates and total amount of conservation laws stay fixed. In BIOMD738, we know that at least Hepcidin is 2-homeostatic (and independent) despite the timeout.

Table 9. Models with more than 12 species tested for ACR, the second number indicate the number of species that verify 2-homeostasis (ACR included). The minimum and maximum value of each total amount of conservation laws is the value computed from the SBML file, divided by 10 and multiplied by 10, respectively. Kinetic rates and volume compartments stay fixed. In BIOMD489, LPS:LBP:CD14: TLR4:TIRAP:MyD88:IRAK4, IkBb_mRNA, IkBe_mRNA, LPS:LBP:CD14:TLR4:RIP1:TRAM:TRIF:TBK/IKKe are detected ACR despite the timeout.

Table 10. Models with more than 10 species tested for homeostasis w.r.t. kinetics rates. The minimum and maximum value of each kinetic rate is given by the SBML file divided by 100 and multiplied by 100, respectively. In BIOMD048, EGF is detected 2-homeostatic despite the timeout. In BIOMD093, SHP2 is detected homeostatic despite the crash.