Identifying parameter regions for multistationarity

Abstract

Mathematical modelling has become an established tool for studying the dynamics of biological systems. Current applications range from building models that reproduce quantitative data to identifying systems with predefined qualitative features, such as switching behaviour, bistability or oscillations. Mathematically, the latter question amounts to identifying parameter values associated with a given qualitative feature. We introduce a procedure to partition the parameter space of a parameterized system of ordinary differential equations into regions for which the system has a unique or multiple equilibria. The procedure is based on the computation of the Brouwer degree, and it creates a multivariate polynomial with parameter depending coefficients. The signs of the coefficients determine parameter regions with and without multistationarity. A particular strength of the procedure is the avoidance of numerical analysis and parameter sampling. The procedure consists of a number of steps. Each of these steps might be addressed algorithmically using various computer programs and available software, or manually. We demonstrate our procedure on several models of gene transcription and cell signalling, and show that in many cases we obtain a complete partitioning of the parameter space with respect to multistationarity.

Fig 1. Running example.

Example network with two species, X₁ and X₂, and three reactions with mass-action kinetics.

Fig 2. Step 1-3 check the assumptions of Corollary 1 and Corollary 2.

In step 4 the function φ_c(x) is constructed and the determinant of M(x) is found. Step 5 is the sign analysis of the polynomial det(M(x)) for $x\in {\mathbb{R}}_{>0}^2$. Step 6 establishes a positive parameterization and finds the polynomial $a\left(\widehat{x}\right)$. Step 7 is similar to step 5, but for $a\left(\widehat{x}\right)$.

Fig 3. Decision diagram of the algorithm.

At each step either the condition is fulfilled or the algorithm terminates indecisively. If that is the case, the corresponding condition might still be verified manually and the algorithm resumed from the next following step.

Fig 4. Two examples describing a hybrid histidine kinase (row 1) and a gene transcription network (row 2).

Column 1: the reaction network; Column 2: the function $a\left(\widehat{x}\right)$ where monomials with coefficients of constant sign (−1)^s are in blue, and those that can have sign (−1)^s+1 are in red; Column 3: parameter conditions for multistationarity; Column 4: Newton polytope where each point corresponds to the exponent vector of a monomial of the numerator of $a\left(\widehat{x}\right)$ (e.g. (1, 2) is the exponent vector of the monomial $x_4{x}_5^2$), blue points are the vertices of the Newton polytope and red numbers indicate the exponents of the red monomials in column 2.