Thalamocortical bistable switch as a theoretical model of fibromyalgia pathogenesis inferred from a literature survey

Fibromyalgia (FM) is an unsolved central pain processing disturbance. We aim to provide a unifying model for FM pathogenesis based on a loop network involving thalamocortical regions, i.e., the ventroposterior lateral thalamus (VPL), the somatosensory cortex (SC), and the thalamic reticular nucleus (TRN). The dynamics of the loop have been described by three differential equations having neuron mean firing rates as variables and containing Hill functions to model mutual interactions among the loop elements. A computational analysis conducted with MATLAB has shown a transition from monostability to bistability of the loop behavior for a weakening of GABAergic transmission between TRN and VPL. This involves the appearance of a high-firing-rate steady state, which becomes dominant and is assumed to represent pathogenic pain processing giving rise to chronic pain. Our model is consistent with a bulk of literature evidence, such as neuroimaging and pharmacological data collected on FM patients, and with correlations between FM and immunoendocrine conditions, such as stress, perimenopause, chronic inflammation, obesity, and chronic dizziness. The model suggests that critical targets for FM treatment are to be found among immunoendocrine pathways leading to GABA/glutamate imbalance having an impact on the thalamocortical system. Supplementary Information The online version contains supplementary material available at 10.1007/s10827-022-00826-8.

where the variables S, T and V represent the mean firing rates of neuron populations belonging to SC, TRN, and VPL, respectively (with the neuron baseline firing rate set to 0), while f (·) is a generic increasing Hill function, g(·) is a generic decreasing Hill function and h(·) is a generic increasing Hill function plus a constant. The strength of the inhibitory interaction (giving rise to both negative feedback loops) is represented by parameter a ≥ 0, which represents the efficacy of the inhibitory GABAergic activity.
Here, we prove that the following qualitative properties hold.
• The system always admits the zero equilibrium (corresponding to baseline activity for all neuron populations), which is always stable.
• When the inhibitory GABAergic activity is absent, i.e. a = 0, the system is cooperative and it may exhibit multistability. The number of equilibria is then odd (including the zero equilibrium).
• For any possible strength of the GABAergic inhibitory interaction a ≥ 0, the equilibria can always be ordered, namely, at each equilibrium the values of all variables (S, T , V ) are larger than the corresponding values at the previous equilibrium. For instance, in the case of three equilibria, the equilibrium values for the variables S, T , and V are ordered as follows: • Which type of behavior can we expect for a small? If for a = 0 there are three equilibria, then the system is bistable and it remains bistable for a small enough. In particular, there are two stable equilibria, the "low" zero equilibrium (S 0 ,T 0 ,V 0 ) = (0, 0, 0) and the "high" equilibrium (S 2 ,T 2 ,V 2 ), while the equilibrium (S 1 ,T 1 ,V 1 ) in the middle is unstable. There exists a threshold a * such that the bistability property remains valid for all 0 ≤ a ≤ a * .
• Which type of behavior can we expect for a large? There is a thresholdā such that, for all a >â, no positive equilibria can exist: the only admissible equilibrium is the zero equilibrium corresponding to baseline neuron activity, (S 0 ,T 0 ,V 0 ) = (0, 0, 0), which is stable.
• Corresponding to the critical value a * , a single positive equilibrium exists, besides the equilibrium at zero, which is shown to be unstable.
An important conclusion can be drawn based on our model: the transition to a pathological state is not continuous, because a threshold a * is expected to lead to an abrupt change in the qualitative behaviour. Precisely, a * corresponds to the transition value above which only the zero (healthy) equilibrium is possible and below which also pathological equilibria exist. For values of a approaching 0, the stable pathological equilibrium becomes the highest possible in terms of the values ofS,T andV .

Model Formulation and Assumptions
Definition 1 A twice differentiable function f (x), defined for x ≥ 0, is a sigmoidal function if f (0) = 0 and f is increasing (for positive values of the argument) and asymptotically constant: lim x→∞ f (x) = c, with 0 < c < ∞. Moreover, the derivative of f is zero both at 0, f (0) = 0, and at infinity, lim x→∞ f (x) = 0, and it has a single maximum for some positivex.
An example of a sigmoidal function is the Hill-type expression considered in the main paper: with p > 1 and positive real parameters α and β.
is a sigmoidal function, namely, if g(0) = c > 0 and g is decreasing (for positive values of the argument) and asymptotically zero: lim x→∞ g(x) = 0. Moreover, the derivative of g is zero both at 0, g (0) = 0 and at infinity lim x→∞ g (x) = 0, and it has a single minimum for some positivex.
An example of a complementary sigmoidal function is the Hill-type expression considered in the main paper: g(x) = γ δ p + x p with p > 1 and positive real parameters γ and δ. Then, the model considered in the main paper belongs to the general class where functions f 1 and f 2 are assumed to be sigmoidal (as per Definition 1, see Figure  1 left), while function k(aT, S) is decreasing in the first argument and increasing in the second (see Figure 1 right).
For a fixed value of T = T 0 > 0, functionk(S) . = k(aT 0 , S) is sigmoidal (as per Definition 1, hence it has the same properties stated before). Conversely, for a fixed value of S = S 0 > 0,k(aT ) . = k(aT, S 0 ) is a decreasing function of aT , asymptotically converging to 0: lim aT →∞k (aT ) = lim aT →∞ k(aT, S 0 ) = 0. We also assume that An example of such function k is the one considered in the main paper: where g is a complementary sigmoidal function (as per Definition 2), and

Model analysis
We are interested in a parametric investigation with respect to a ≥ 0 (in the main paper, we focus on 0 ≤ a ≤ 1).

Boundedness
Since the functions f 1 , f 2 and k are bounded (they cannot exceed a maximum value, regardless of the values of their arguments), the solutions of the system are bounded due to the presence of the linear terms.

Equilibria
Let us consider the set of possible equilibria. The zero vector (S 0 ,T 0 ,V 0 ) = (0, 0, 0) is always an equilibrium, associated with the baseline neuron activity (healthy state). This can be immediately checked by substituting the zero vector in the system equations (4)-(6): then, the resulting derivatives become zero. We wonder whether there are other equilibrium points. The equilibria are derived by settingṠ = 0,Ṫ = 0,V = 0 and finding the solution of the resulting system of algebraic equations: As anticipated, (S 0 ,T 0 ,V 0 ) = (0, 0, 0) is always a solution. In general, from (9) we get where we consider V as a variable and a as a parameter. The roots of function Ψ in V correspond to the possible equilibrium values, for a given a.  Before proceeding with our analysis, we introduce the notion of critical root.
Definition 3 Given the equation f (V ) = 0, a root V * is critical if, besides being a zero of the function, f (V * ) = 0, it is a zero also of the function derivative: f (V * ) = 0. For our considered system, we call critical equilibrium an equilibrium corresponding to a critical root of function Ψ(a, V ).
We may have a critical root V * of Ψ(a, V ) for some critical value of a, a = a * . In fact, we can have the following possible cases, depicted in Fig. 2.
• For large values of a, there is only the root at zero (cf. the red curve in Fig. 2).
• For small values of a, in addition to the root at zero, there are multiple roots (cf. the blue curve in Fig. 2), all included between the extremal positive roots of Ψ(0, V ), corresponding to a = 0 (cf. the black curve in Fig. 2, with extremal positive rootsV andV ).
• There is a critical intermediate value, a = a * , for which a critical root exists (cf. the green curve in Fig. 2), in addition to the root at zero.
We can formalize these findings in the following propositions.
Proposition 1 Assume that there are no critical roots. Then, the number of nonnegative roots of equation Ψ(a, V ) = 0 is odd.
because both f 1 (0) = 0 and f 2 (0) = 0. Hence, for any fixed a, function Ψ(a, V ) is negative in a right neighborhood of 0. For V large, function Ψ(a, V ) tends to −∞, because function k is a composition of bounded functions, while function (−V ) tends to −∞. Therefore, the curve can intersect the abscissa axis an odd number of times overall (counting also the intersection at 0). If, for some value a = a * , there is a critical root V = V * (see the green curve shown in Fig. 2), Proposition 1 still holds provided that the roots are "counted with their multiplicity". Since Ψ(a * , V * ) = 0, d dV Ψ(a * , V * ) = 0 and the second derivative is non-zero, the Taylor expansion is: . . , and the critical root has multiplicity 2; the roots have a multiplicity corresponding to the index of the first non-zero derivative. Later on, we will analyse the stability properties of the critical equilibrium, corresponding to a critical root, and conclude its instability. In case no positive equilibria exists for a = 0, namelyV =T =S = 0 is the unique equilibrium when a = 0, then no positive equilibrium can exist also for a > 0.
Proof: The result follows immediately from the fact that, being function k decreasing in the first argument, increasing the value of a shifts the curve down-word. Therefore, since Ψ(a, V ) ≤ Ψ(0, V ) and Ψ(0, V ) < 0 for V >V , the rootsV a of Ψ(a, V ) = 0 are smaller thanV :V a <V . In turn,S a = f 1 (V a ) < f 1 (V ) =S andT a = f 2 (V a +S a ) < f 2 (V +S) = T , because both f 1 and f 2 are increasing functions. The proof for the lower bounds is identical.

Stability analysis
We consider the following definition of stability.
Definition 4 An equilibrium is stable if the linearized system around that equilibrium has eigenvalues with strictly negative real part. Otherwise, the equilibrium is unstable.
In principle, this definition is restrictive, since there are examples of systems whose linearization has eigenvalues with non-positive real part, including eigenvalues with zero real part, and whose trajectories converge to the equilibrium if the initial condition is close enough. However, we do not contemplate this fragile situation as a stable equilibrium, hence we rule out eigenvalues with zero real part, because infinitesimal perturbations can make their real part positive.
We now prove that the healthy equilibrium at zero, corresponding to baseline activity for all neuron populations, is always stable.
Proof: Since the derivatives of functions f 1 , f 2 and k are zero when computed at the origin, the linearized system at the origin is associated with a diagonal state matrix whose eigenvalues are the diagonal elements, which are real and negative. Hence, the equilibrium is stable.
The theory that follows allows for any (odd) number of equilibria. However, to avoid cumbersome notations, we assume that the equilibria can be at most three: the possible roots of ψ(a, V ), for a given a, are of which only the zero root surely exists. Correspondingly, the values for S and T are This proves the ordering among the equilibria. We can also prove the following result.
Proposition 4 For a = 0, if there are three ordered equilibria, the zero equilibrium (S 0 ,T 0 ,V 0 ) = (0, 0, 0) and the "high" equilibrium (S 2 ,T 2 ,V 2 ) are stable, while the medium equilibrium (S 1 ,T 1 ,V 1 ) is unstable. There exists a threshold a * > 0 such that this pattern is preserved for all a < a * .
Proof: The result immediately follows from the fact that, for a = 0, the system is monotone (see [1] and [2]) and therefore enjoys special properties: its equilibria are ordered and alternating, stable and unstable. The last claim follows by continuity.

Equilibria for large a
Finally, the next proposition states that, for a large enough, there are no positive roots for Ψ(a, V ) = 0.
Proposition 5 There exists a thresholdâ such that, for a >â, the only possible equilibrium is the one at zero.
Proof: Assume that two positive roots exists for Ψ(0, V ) = 0, denoted asV 0 andV 0 , with 0 <V 0 <V 0 . For any a, all other positive roots must be in the interval [V 0 ,V 0 ] (whose extrema are the black dots shown in Figure 2), according to Proposition 2. Now, by contradiction, assume that no such thresholdâ exists. Hence there exists a sequence a j → ∞ for which a positive rootV a j exists, which must necessarily be within the interval [V 0 ,V 0 ]. The sequence of such roots {V a j } has an accumulation pointV * . Without loss of generality, we can assume that {V a j } converge toV * (if not, we just need to consider a proper converging sub-sequence). Then, since Ψ(a j , V a j ) = 0 and V a j →V * , necessarily Ψ(a j , V * ) → 0 for some j. On the other hand, by assumption k(a j T, S) converges to 0 as a j → ∞. Since , then V * must be 0, which is not possible because V * ≥V 0 > 0. The result is therefore proven by contradiction.

Critical equilibria
When a critical equilibrium is present, then we can show that it is unstable.
Proposition 6 Assume that, besides the equilibrium at zero, there exists a critical equilibrium. Then, the critical equilibrium is unstable.
Proof: Let us rewrite the derivative of Ψ(a, V ) as follows: The derivative is zero by assumption, because we are considering a critical equilibrium with S * = f 1 (V * ) and T * = f 2 (V * + S * ). Now consider the Jacobian J of the system factorized as follows: which is in line with the definition of critical equilibrium. This means that the Jacobian is singular, hence it admits λ = 0 as an eigenvalue, which implies the instability of the equilibrium, according to our definition.
Although the scenario with a critical equilibrium is unlikely to occur, because it happens for a specific value of a * , the analysis is interesting also to understand the case in which two equilibria are very close: in this case we could have, in principle, stability of the upper equilibrium, but its domain of attraction would be very small. This implies that, for the upper (pathological) equilibrium to be stable with a large domain of attraction, the value of a must be sufficiently small.