Transition from Thermokinetic to Chemical Instabilities in a Modified Sal’nikov Model

Abstract

Traditionally, thermokinetic and chemical oscillations have been studied independently, but in cellular media recent studies have shown that cell’s temperature is not uniform. Thus, on this context it is possible to inquire about the influence of thermal effects on chemical oscillatory behavior and vice versa. To this end, in this paper we propose a dynamical model based on a modified Sal’nikov oscillator that can address both kinds of oscillatory behavior (thermokinetic and chemical). We found that the system modeled can jump from thermal oscillations to mixed chemical-thermal oscillations through a transition or bifurcation parameter. Thermokinetic oscillations are well defined in a limit cycle, while chemical-thermal oscillations appear in the form of a burst. The model could be useful in explaining biochemical energy recovery under cellular stress conditions.

Keywords:
chemical instabilities; dynamical transitions; thermokinetic instabilities; biochemical oscillations

Introduction

The kinetics of chemical oscillations and thermal oscillations have been extensively studied in an independent manner. Studies regarding the former are based on the well-known Belousov-Zhabotinskii reaction,^1-3 while the latter have been addressed by applying the principles of the Sal’nikov model to combustion reactions.^4-7 However, while this gap is not necessarily essential in the chemical engineering field, since temperature can be controlled in several chemical processes,⁸ in the biological context of a living cell it could be relevant at the moment of examining the biochemical reactions that take place, first, in non-isothermal conditions and second, in the presence of chemical feedback. The first condition is known due to recent advances in methods intended to map intracellular temperature, which have made it possible to measure temperature changes in individual cells.^9-11 In fact, Takei et al.¹² used a fluorescent thermometer to demonstrate that Ca²⁺ oscillations inside HeLa cells correlated well with temperature increments. Additionally, Arai et al.¹³ reported increments in the intracellular stem temperature caused by Ca²⁺-adenosine triphosphate-ase (ATPase) pumping. The second condition results from chemical feedback inherent to the biochemical cycles, for example in the form of chemical autocatalysis.¹⁴ Thus, it is feasible for variable non-isothermal conditions and chemical feedback to occur at the same time, and the system behavior is expected to be influenced by both sources of instability.

Initially, thermokinetic oscillations were of primary interest in the chemical engineering field, where the role of temperature in reactor operation is fundamental, especially if exothermic chemical reactions are involved.^15,16 These oscillations are particularly relevant in cases where uncontrolled temperature increases represent a security risk, since thermal runaway development could be a lethal hazard.¹⁷ Additionally, several examples of thermal oscillatory behavior can be found in combustion reactions, traditionally modeled as the Sal’nikov thermokinetic oscillator.^4-7 According to the mechanism of the model, the Sal’nikov oscillator is driven only by temperature changes, and is said to be a pure thermal oscillator.

In addition to pure thermal oscillators, sometimes systems exhibiting only chemical instabilities are prone to displaying thermokinetic oscillations if they reach an operation point over non-isothermal conditions, either because the chemical reaction itself generates heat or because it is subjected to an external heating bath. In both cases, heat flux is incorporated through the energy balance equation with nonlinear terms in the Arrhenius rate constants.^4,6 However, it is common to regard the thermal effect as secondary when chemical instabilities are present, and to introduce a thermostatic control that guarantees an average isothermal reaction. This approach could hide interesting features exhibited by chemical oscillators, for instance when analyzing the Belousov-Zhabotinskii reaction, knowing in advance that some reaction steps are highly exothermic.^18-21 This particular scenario, where chemical and thermal instabilities are present, has led to questions about the source of the oscillatory behavior and how chemical contributions can be distinguished from thermal ones.²²

A model able to reproduce both chemical and thermal oscillations could be valuable in a biological context. To mention one hypothetical case, muscle contraction under normal conditions requires ATP provided by oxidative metabolism, but at the onset of a high-intensity exercise is exceeded its capacity to provide energy. Stored phosphocreatine cannot sustain such continuous intense activity either, giving way to the glycolytic pathway as a main energy source.²³ But even glycolysis is not able to bypass muscle fatigue. When muscle is exhausted, it is necessary to recover the energetic state, and, in fact, muscle fatigue is a mechanism to protect the cell from low energy levels.²⁴ With the purpose of resynthesize ATP, biochemical instabilities characteristic of glycolytic pathway¹⁴ and non-isothermal medium where muscle glycolysis take place could be coupled. The coupling between thermal gradient and glycolytic oscillations to force the Ca²⁺-ATPases (associated with a heat-generating process)^25,26 to synthesize ATP could be an emergency mechanism to address cellular stress (see Figure 1).

Figure 1
Hypothetical biological context where chemical oscillations from glycolysis could coincide with thermal gradients arising as a consequence of muscle exercise. At this point, heat entering Ca²⁺-ATPase could force this enzyme to resynthesize ATP with the ADP coming from phosphofructokinase (PFK) step of glycolysis.

In this paper, we propose a dynamical model based on a modified Sal’nikov oscillator that can address both kinds of oscillatory behavior (thermokinetic and chemical). Numerical stability analysis confirms that in a definite region of interest, through a bifurcation parameter γ, the model can jump from thermokinetic instabilities to chemical-thermokinetic instabilities. Rate constants of the chemical network are considered of Arrhenius type. But, instead of forcing temperature changes with an auxiliary mathematical function, it is included an additional energy balance equation. This equation depends both on the exothermicity of the chemical reactions as well as the thermal gradient between reaction medium and the surroundings.

Methodology

The ultimate purpose of this paper is to demonstrate a possible alternation between chemical and thermokinetic oscillations in a modified Sal’nikov model. From this perspective, it is necessary to modify progressively the original two-step Sal’nikov model to add an intermediate step (three-step model), and finally a chemical autocatalysis (unified model). The alternation between original model and the model with chemical feedback is mediated by a bifurcation parameter γ.

Sal’nikov was the first to propose a simple model that could generate temperature oscillations in non-isothermal conditions.⁴ Although the model was intended to explain the occurrence of cool flames observed during the oxidation of hydrocarbons, it has been continuously studied, since it involves thermal instability as a source of kinetic oscillations, which coincides completely with chemical and energy principles.²⁷

The Sal’nikov scheme consists of two consecutive first-order reactions followed by an exothermic reaction. It closely follows the model formulated by Gray et al.⁷ in the context of thermal combustion, with a precursor P generating reactive species Y which decomposes through an exothermic reaction to produce inert product B. It is assumed that all the heat output (enthalpy of the reaction, Q) is associated with the second reaction. This set of reactions is represented as:

(1) $P\stackrel{k_1,E_1}{\to }Y\left(Q_1=0\right)$

(2) $Y\stackrel{k_2,E_2}{\to }\left(Q_2\ne 0\right)$

where the temperature dependence of the two reactions takes Arrhenius form, k₁ and k₂ are the rate constants, and E₁ and E₂ are the energies of activation. Step 1 (equation 1) is taken to be thermoneutral and to have a rate that is independent of temperature (i.e., zero activation energy, E₁= 0); it is also assumed that the precursor reactant undergoes a first-order decay. We refer to this model as Sal’nikov’s two-step model in order to distinguish it from the modified Sal’nikov three-step model, or unified model (in addition to the Y intermediate, this model includes a second intermediate, Z), which will be introduced later.

The set of reaction rates ν₁ and ν₂ associated with equations 1 and 2 are written as:

(3) $v_1=p_0{k}_1\exp \left(-k_{1}t\right)$

(4) $v_2=k_{2}y$

where p₀ is the initial reactant concentration, t is time and y is the concentration of reactive species.

Reaction rate ν₁ is expressed as the exact solution for the dynamics of equation 1, and is presented in this form as the most general case. The set of governing equations will be:

(5) $\frac{\mathrm{dy}}{\mathrm{dt}}=v_1-v_2$

(6) ${\mathrm{nC}}_p\frac{\mathrm{dT}}{\mathrm{dt}}=Q_2{v}_2-S\chi \left(T-T_a\right)$

In this model, rate constant k₂ takes the form:

(7) $k_2\left(T\right)=k_2\left(T_a\right)\exp \left(\left(E_2/R\right)\times \Delta T\right)$

where $\Delta T=1/T_a-1/T$, because thermal dependence is linked to the conversion $Y\stackrel{k_2}{\to }B$. C_p is the heat capacity of the mixture, n is the total number of moles, S is the reactor surface area, χ is the surface heat transfer coefficient and T_a is the surrounding (environmental) temperature. Heat loss is modeled as a simple Newtonian heat transfer, assuming a well-stirred solution. Table 1 summarizes the parameter values used to simulate system dynamics, according to equations 1 and 2.

Table 1
Representative values of Sal'nikov physicochemical quantities

Physicochemical quantity	Symbol	Value
Volume / m3	V	0.001 (cylindrical reactor with h = 0.15 m and r = 0.046 m)
Surface area / m2	S	0.05
Surrounding temperature / K	Ta	400
Activation energy / (kJ mol-1)	E	166
Heat capacity / (kJ K-1 mol-1)	Cp	0.029
Surface heat transfer coefficient / (kW m-2 K-1)	c	0.055
Change in the enthalpy of exothermic step / (kJ mol-1)	Q	400
Initial reactant concentration / (mol m-3)	p0	1
Rate constants	k1	0.01 s-1
	k2a	400 M-1 s-1
Gas constant / (kJ mol-1 K-1)	R	0.008314

• ^a
  Constant k₂ could obey to a first order or second order kinetics. Although, its value is assumed to be the same.

The different models were characterized with a classical linear dynamic analysis focused on verifying oscillatory behavior for a set of values of the interest parameters (k_2,0, k₂ and χ). Rate constant k_2,0 sub-indice indicate that it is the value of the rate constant at the surrounding reference temperature or, in other words, k₂(T_a). For the the three-step model, k_3,0, the third step rate constant, will be taking into account to represent stability analysis results. The previous procedure was implemented in Matlab^® R2014b²⁸ in order to construct a representative qualitative plot (2D or 3D, depending on the model’s dimensions) that would indicate the stability of the system. Matlab’s script gave a set of points that satisfied the condition of an imaginary part different from zero in the eigenvalues entries of the stability matrix, related here with unstable states, and represented in an appropriate space defined by the interest parameters.

The dynamics of Sal’nikov’s two step model has been well-established.^7,29,30 Starting with a set of known values for the parameters that guarantee thermokinetic oscillations, the addition of an intermediate step should not modify neither stability analysis results nor oscillatory behavior. Then, the autocatalytic step is fully activated and the associated stability surface is overlapped with the correspondig to the two-step model to find a common region where without changing parameters values is yet possible observe both chemical and thermokinetic oscillations.

Results

Sal’nikov’s three-step model

An additional intermediate step between steps 1 and 2 (equations 1 and 2, respectively) is inserted in order to preserve the original schema while gaining the possibility of later manipulating the new step to add the chemical feedback. Before proceeding with this subtle change, it is necessary to demonstrate that the original Sal’nikov model remains the same when inserting a third intermediate species, Z. This condition is satisfied if the stability surface obtained from the modified Sal’nikov model is consistent with the 2D diagram of the original model (see Figure 2b).

Figure 2
Oscillatory regions for two-step or three-step model alongside k₂, k_2,0 or k_3,0 (k₂ is treated as a variable), and χ, when other values are those found in Table 1, under the conditions Q₁ = Q₂ = 0 and Q₃ ≠ 0. (a) The black dots delimit the oscillatory region in the space χ - k_2,0 if only the two-step model is considered (k_2,0 = k₂(T_a)); (b) the surface is actually the projection of this area extended over values of k₂ (three-step model) and becomes a frontier between stable behavior (over) or unstable behavior (underneath). The region starting from 0.04 kW m^-2 K^-1 is only shown inasmuch as it takes the form of a plane, until 0.2 kW m^-2 K^-1, where the surface disappears because all points drive the system to a stable state (k_3,0= k₃(T_a)).

After adding the intermediate species Z, the set of equations 1-2 becomes 8-10, as shown below:

(8) $P\stackrel{k_1,E_1}{\to }Y\left(Q_1=0\right)$

(9) $Y\stackrel{k_2,E_2}{\to }Z\left(Q_2\ne 0\right)$

(10) $Z\stackrel{k_3,E_3}{\to }B\left(Q_3\ne 0\right)$

The same considerations stated for the original model are preserved, though due to the presence of the new step it is possible to explore the model’s response when taking the exothermic reaction as the second one (Q₁= Q₃ = 0 and Q₂ ≠ 0) or the last one (Q₁ = Q₂ = 0 and Q₃ ≠ 0). Stability analysis reveals the expected agreement between the two-step and three-step models, since the 3D surface is the projection of the stability diagram over k₂ values in the 2D space defined by χ - k_2,0, when it is assumed that Q₁= Q₂= 0 and Q₃ ≠ 0 (see Figure 2b).

Once it has been demonstrated that the two-step and three-step models are essentially the same in terms of their stability analysis, it is possible to study the response of the three-step model when heat output is placed completely on the second reaction (equation 6) or third reaction (equation 7). According to Figure 3a or 3b, depending on the step selected to be the most exothermic, there is a region where each particular species involved will oscillate (Y and Z simultaneously when Q₂ ≠ 0, or only Z when Q₃ ≠ 0). This is relevant in a context in which there is an interest in having one oscillatory behavior at the output of the system while the other intermediates remain stable.

Figure 3
Effect of assuming heat output in steps 6 or 7 on Sal’nikov’s three-step model. Upper curves of each figure correspond to Z (continuous line) and Y (dashed line) species and lower plot corresponds to temperature changes. (a) For k₁ = 0.01 s^-1, k_2,0 = 0.8 s^-1 and k₃ = 100 s^-1, the Y and Z intermediates oscillate if Q₁ = Q₃ = 0 and Q₂ ≠ 0; (b) only the Z intermediate oscillates if Q₁ = Q₂ = 0 and Q₃ ≠ 0, when k₁ = 0.01 s^-1, k₂ = 1 s^-1, and k_3,0 = 0.2 s^-1. The values not mentioned here are the same as those found in Table 1.

Sal’nikov’s unified model: combining chemical and thermokinetic oscillations

The Sal’nikov models, presented up to this point, illustrate only the thermokinetic aspect of the phenomena in such a way that it is not possible to establish the influence of chemical instability, since there is no chemical feedback loop. These independent portrayals should be unified in order to get the full picture. This can be achieved by considering a simpler Sal’nikov model that incorporates both cases in a continuous manner through a transition parameter, making it feasible to distinguish the transition from chemical to thermal oscillations and vice versa, and to determine under what conditions this bifurcation parameter triggers transitions.

Again, we consider the simple global reaction P → B, which for this case will be the sum of three reaction steps (equations 5-7). This is no different from the previous model. However, adding a “switch” parameter, γ, opens the possibility to gain an autocatalytic step, if γ = 1, or permits the recovery of the three-step model when γ = 0. The model is based on equations 8-10:

(11) $P\stackrel{k_1,E_1}{\to }Y$

(12) $Y+\gamma Z\stackrel{k_2,E_2}{\to }\left(1+\gamma \right)Z$

(13) $Z\stackrel{k_3,E_3}{\to }B$

Following auxiliary equations for reaction rates, it is assumed that:

(14) $v_1=p_0{k}_1\exp \left(-k_{1}t\right)$

(15) $v_2=k_2{\mathrm{yz}}^{\gamma }$

(16) $v_3=k_{3}z$

In this case, thermal dependence is linked to the conversion $Z\stackrel{k_3}{\to }B$, where k₃ takes the form:

(17) $k_3\left(T\right)=k_3\left(T_a\right)\exp \left(E_3/R\times \Delta T\right)$

Finally, the set of governing equations is:

(18) $\frac{\mathrm{dy}}{\mathrm{dt}}=v_1-v_2$

(19) $\frac{\mathrm{dz}}{\mathrm{dt}}=v_2-v_3$

(20) ${\mathrm{nC}}_p\frac{\mathrm{dT}}{\mathrm{dt}}=Q_3{v}_3-S\chi \left(T-T_a\right)$

In contrast with the stability analysis of the previous model, this model exhibits a 3D surface with two unstable regions (see Figure 4, γ = 1 surface). One is defined for small values of χ (0.01-0.2 kW m^-2 K^-1) and the other for values of χ higher than 0.2 kW m^-2 K^-1. In the first range (see Figure 5, plot at plane χ = 0.01), oscillations appear in the form of a long burst both in the Z concentration and the temperature. In the second region (see Figure 5, plot at plane χ = 0.5), the parameters force the system to exhibit an initial peak in the concentration of the species involved, as well as temperature, finally evolving to a fixed value. This surface summarizes dynamical behavior when the case γ = 1 is considered, i.e., oscillations due to the presence of a chemical feedback loop.

Figure 4
Stability surface plots for Sal’nikov’s unified model alongside k₂, k_3,0 (explored until k_3,0= k₃(T_a) = 3 s^-1) and χ, when other values are those found in Table 1. γ = 0 surface was previously plotted in Figure 2a. Each surface represents the frontier of chemical instabilities: all points above them drive the system to global qualitative stable states and the points underneath represent unstable states (rate constant k₂ has units of s^-1 if γ = 0 or M^-1 s^-1 if γ = 1). It is of special interest the common region where both γ = 0 and γ = 1 share the same set of unstable state points.

Figure 5
Transition from thermokinetic and chemical oscillations (χ = 0.01 kW m^-2 K^-1) to pure chemical oscillations (χ = 3 kW m^-2 K^-1) when k₁ = 0.01 s^-1, k₂ = 400 M^-1 s^-1 and k_3,0 = 0.7 s^-1 for the unified model (γ = 1). Upper curves of each plane correspond to Z (green) and Y (blue) species and lower plot correspond to temperature changes. Time scale is the same for a particular plane, and for this reason is only represented once. Values not listed here are those found in Table 1.

In this way, we have shown the thermokinetic and chemical effects on a modified Sal’nikov model, though these are still partial understandings of a generalized case. A complete picture of the system dynamics is provided by superimposing the stability surfaces corresponding to the values of γ = 0 and γ = 1. Figure 4 combines the surface for γ = 0 (three-step model, equations 5-7) with the surface for γ = 1, explained above. It is clear that for small values of χ both surfaces share the same stability space where it is thus possible to obtain oscillations of thermokinetic nature (γ = 0, without chemical feedback) or oscillations in the presence of a chemical feedback loop (γ = 1, with thermokinetic influence). Thus, by means of the stoichiometric bifurcation parameter γ, the system jumps from one oscillatory condition to the other, as shown in Figure 6 (keeping all parameters the same, except for γ, it is possible to get sustained oscillations (Figure 6a) or oscillatory bursts (Figure 6b)).

Figure 6
Comparison of system dynamics for χ = 0.05 kW m^-2 K^-1 (the χ value is picked from the common region in the surfaces from Figure 4); k₁ = 0.01 s^-1, k₂ = 400 M^-1 s^-1 and k_3,0 = 0.7 s^-1. Values not listed here are those found in Table 1. Upper plots correspond to Z (green) and Y (blue) species. Lower plots correspond to temperature changes. (a) When γ = 0, the system exhibits sustained oscillations both in the Z species and temperature T; (b) for the same set of parameter values, and using only γ = 1, sustained oscillations are substituted for a long burst in the temperature and the Z chemical concentration.

Even when Figure 4 supports the transition from pure thermokinetic oscillations to mixed oscillatory behavior, it is necessary to study the predominance either of thermal effects or chemical feedback loop. According to Figure 5 (which illustrates changes in the variables of interest for a combination of the parameters taken from the common space of surfaces γ = 0 and γ = 1), if the heat transfer coefficient chosen is in the range of values between 0.01 and 3 kW m^-2 K^-1, the unified system output still oscillates, either due to the thermal effect if χ is low or due to chemical feedback if χ is high. Furthermore, if thermal effect is removed completely from equations 11-13, i.e., simulations are performed without equation 13, the system will be reduced to the isothermal case, reflecting only chemical autocatalysis. As shown in Figure 7, the isothermal curve is similar in shape to the non-isothermal (Figure 5, plane at χ = 3 kW m^-2 K^-1), and difference is because in the non-isothermal case yet remains a heat contribution from reaction’s exothermicity.

Figure 7
Comparison between isothermal and non-isothermal Sal’nikov system. For the isothermal case, only it is assigned a high value for χ in the equations 11-13. For the non-isothermal case, the system is converted in a pure chemical oscillator by removing the thermal feedback from equation 13. The shape of non-isothermal oscillations follows the dynamic of the pure chemical oscillator. Upper curves of each figure correspond to Z (continuous line) and Y (dashed line) species and lower plot corresponds to temperature changes.

The main difference between Sal’nikov’s two-step model and the unified model consists in the addition of a chemical feedback loop, where the transition is mediated by the γ value. This parameter appears in a natural way because a stoichiometric factor presents itself for the new chemical species, Z. When considering a particular mechanism reaction, such as the one stated in equations 8-10, γ should only take the discrete values 0 or 1. This stoichiometric balance remains for a single reaction chain, but when a complex reaction chain is considered, for instance those occurring inside a living cell, it is perfectly plausible to expect a continuous set of γ values ranging from 0 to 1. Figure 6 only shows the system’s output when considering the cases of γ = 0 or γ = 1, and Figure 8 depicts how the unified model varies for different γ values when the model’s input is fed with particular parameters taken from Figure 5. This type of oscillation thus depends on the stoichiometric bifurcation parameter γ, and the unified model addresses changes associated with different γ values.

Figure 8
Expected output from unified model after varying the stoichiometric factor (bifurcation parameter) γ between 0 and 1 (parameter set values are the same as those considered in Figure 7). Upper curves of each plane correspond to Z (green) and Y (blue) species and lower plot correspond to temperature changes. Time scale is the same for a particular plane, and for this reason is only represented once. Assigning different values to γ confirms the fact that the unified model is not only limited to a single reaction chain occurring inside a hypothetical reactor, but could be extended to a biochemical context, where γ can take different values according to the evolution of the biochemical reaction network over time.

In order to characterize the stability of the model for different γ values, in the Figure 9 phase plane plots are constructed between temperature and Z concentration. On this way, closed trajectories are found for all γ considered. Cycle limits are well defined for γ < 0.8, but for γ near to 1, though initially the system enters in a closed loop, gradually gives way to damped oscillations. Additionally, near to γ = 1, amplitude of the Z oscillations tend to increase, while concentration of Y becomes not self-regulated.

Figure 9
Phase plane plots between temperature and Z concentration (parameter set values are the same as those considered in Figure 7). For the tested values of γ, the system enters in a cycle limit trajectory, except for γ = 1, where system’s dynamics corresponds to damped oscillations. Trajectories tend to increase in amplitude near to γ = 1.

Discussion

It has been known for years that thermal oscillations occur in chemical reactions, with autocatalysis as a source of instability, in both isothermal reactors and adiabatic reactors,^16,17,31 the same type of oscillations occur in chemical processes that are known to be thermokinetically unstable.^4-7 Both kinds of oscillatory behavior (thermal and chemical) have been shown in this paper using a feasible dynamic model, which we have called the Sal’nikov unified model (see equations 8-10). Figure 8 summarizes a representative set of plots obtained from the numerical results, which illustrates oscillations in temperature and oscillations in chemical concentrations. The unified model was built from the well-known Sal’nikov model, by adding an advantageous intermediate reaction step (a third step, since the original model only has two) to force variable stoichiometry by means of a bifurcation parameter (Z species stoichiometric factor, γ). There are two noteworthy features of this model: first, different values for the (bifurcation) transition parameter γ drive the system to exhibit thermal oscillations (γ = 0) or chemical and thermal oscillations (γ = 1), without modifying kinetic rate constants or reactor parameters; second, all magnitudes associated with the system’s parameters.^32,33

In order to distinguish thermokinetic oscillations from chemical oscillations, and verify that effectively the model can display chemical instabilities, the heat input through χ was restricted from equation 13. After increasing the value of χ, the system’s ability to retain external heat decreases and external temperature variations no longer affect the system’s temperature. This particular condition evinces chemical damped oscillations. A further comparison with the non-isothermal case, where thermal instability was removed, i.e., running simulations without equation 13, confirmed that non-isothermal oscillations at χ = 3 kW m^-2K^-1 are mainly driven by chemical instabilities (see Figure 5). Difference in amplitude is due to the additional contribution from exothermicity of the $Z\stackrel{k_3}{\to }B$ step (Q₃ν₃ term in equation 13). χ represents itself another bifurcation parameter, but here is merely changed for studying the relative contribution of chemical and thermokinetic instabilities. Fixing all the relevant parameters to the common space where both γ = 0 and γ = 1 surfaces share the same points that potentially could drive the system to display an oscillatory behavior, guarantees that by only changing γ the system will jump from oscillations with chemical origin to oscillations with thermokinetic origin.

Additionally, transition parameter γ can be understood physically if we place our system in a complex reaction context, where other processes that dynamically involve the Z species occur, or in other words, processes that can benefit from a temperature-dependent step in the presence of a chemical feedback loop. This could be the case of skeletal muscle contraction response considered under stressing conditions or muscle fatigue (see the review by Allen et al.).³⁴

From such cellular scenario emerges two hypothetical candidate processes that can benefit one from each other. On one hand, a biochemical process providing energy, the glycolytic pathway, and capable of generating oscillations; and a non-isothermal media, due to the continuous heat contribution by muscle sarcoplasmic Ca²⁺-ATPases (fibers which are capable of faster contractile cycle tend to have a higher density of Ca²⁺-ATPases).³⁵ A closer look on both oscillatory glycolysis and Ca²⁺-ATPase heat generation would show how they could fit into the scheme of unified model.

Glycolytic oscillations are noted in the list of biochemical oscillators because they have been extensively studied in yeast cells, becoming a model of study. Classically, this type of oscillation has been related to the negative feedback loop exerted by ATP/ADP ratio over the phosphofructokinase enzyme (PFK).^14,36 Higgins³⁶ and Sel’kov³⁷ proposed first simplified models that reproduced glycolytic oscillations observed in yeast cells. The non-dimensional equations that describe oscillations are $\frac{\mathrm{dx}}{\mathrm{dt}}=-x+\mathrm{ay}+x^{2}y$ and $\frac{\mathrm{dy}}{\mathrm{dt}}=b-\mathrm{ay}-x^{2}y$, where x and y are concentrations of ADP and fructose-6-phosphate, and a,b > 0 are kinetic parameters.³⁸ A remarkable feature of the original, non-dimensional Sel’kov’s model, was the presence of a stoichiometric parameter that controls the number of molecules that form a complex with phosphofructokinase (PFK). By means of this factor, glycolysis could exhibit different dynamic behaviors. Additional to intrinsic characteristics of the model, external features, like intracellular temperature changes, have been recognized as a feedback control. For example, Mair et al.³⁹ studied the entrainment of yeast glycolytic oscillations by applying thermal pulses and thermal gradients. They found that temperature acts as a controller of glycolytic oscillations, being the PFK enzyme sensitive to thermal changes. Moreover, when the glycolytic system was placed under a thermal gradient, it was induced nicotinamide adenine dinucleotide (NADH) traveling waves.

In the previous example, it is assumed an Arrhenius-type temperature dependence of the kinetic rate constants involved in the chemical reaction network. It is supposed that when changing the temperature, the rate of kinetic reaction changes. But this is an implicit fact, and the non-linearities in the mass-action kinetics are the responsible for the emergence of oscillations. In our case, this is also true if looking at the system when γ = 1, because equations 11-12 describe a system with quadratic autocatalysis. But, additionally to the implicit Arrhenius dependence, a third equation (equation 13) explicitly dictates the shape of the thermal profile. This is subtly different to force the chemical oscillator only to response to a thermal profile, but, as we propose, the system itself incorporates thermal oscillations. As we show in Figure 4, numerical analysis confirms the existence of a region for the parameters of interest where there are chemical and thermokinetic instabilities.

Even if glycolytic cycle is not forced to show oscillations, or is not oscillating, glycolytic oscillations are prone to appear when energy state must be preserved, being one interesting case the ATP depletion induced by intense exercise in skeletal muscle fibers. It has been demonstrated that glycolysis rate in skeletal muscle is controlled by factors related to energy state.⁴⁰ Additional evidence from isolated rabbit ventricular myocites shows that glycolysis can cause periodic oscillations on ATP levels if oxidative phosphorylation and creatine kinase system cannot buffer the cellular ATP/ADP ratio.⁴¹ This special condition guarantees biochemical oscillations, a half part of the Sal’nikov unified model.

The other half part must be a temperature generating process. In the proposed scenario, evidence of temperature-dependent process corresponds precisely to the thermal phenomena inherent to the Ca²⁺ pumping by sarcoplasmic Ca²⁺-ATPase. At this respect, Arai et al.¹³ recorded the dynamics of production of heat from Ca²⁺-ATPase at endoplasmic reticulum of HeLa cells using a molecular fluorescent probe. A gradient temperature, the distinctive element upon which the Sal’nikov unified model is built, stems from the continuous working of the pump. In the particular case of muscle contraction, the active transport of calcium ions by Ca²⁺-ATPases evinces thermal forces and heat fluxes converging in a biochemical process.^25,26,42

In order to recover its energy state, biochemical oscillations and temperature gradients provide a tentative answer to the regulatory response adopted by a biological system. Thus, instead of expecting ATP production to occur solely from a biochemical pathway tied to chemical instabilities, it could be the result of coupling thermal force and chemical machinery, as stated by non-equilibrium thermodynamics.⁴³ In fact, Lervik et al.⁴⁴ used non-equilibrium molecular dynamics simulation to show that Ca²⁺-ATPase is able to sustain large thermal gradients across its structures.

Here we hypothesize that to contribute with the energetic recovery, ATP could be resynthesized by the Ca²⁺-ATPase pump working backwards,^45,46 i.e., producing ATP instead of using it to uptake calcium (see Figure 1). In this way, thermal gradient is used to support the chemical production of a critical specie. Evidence of ATP synthesis from heat fluxes is discussed by Müller,^47-49 who proposes it as a basic mechanism used by the first living organism to transform energy into a usable chemical form. Although this work does not discuss archaic cells, this kind of primitive mechanism could be the cell’s answer when it needs to restore ATP levels.

Conclusions

In this work, we presented a mathematical model based on the Sal’nikov oscillator. The model is proposed in order to address both kinds of kinetic instabilities, thermokinetic and autocatalytic. For this, we added to the Sal’nikov oscillator a third step of reaction that includes a stoichiometric coefficient as a bifurcation parameter. After adjusting all set of parameters, numerical modeling of this modified Sal’nikov oscillator shows that the time series exhibit transitions from thermal oscillations to mixed chemical-thermal oscillations when the stoichiometric coefficient varies between 0 and 1. Thermokinetic oscillations are well defined in a limit cycle, while chemical-thermal oscillations appear in the form of a burst. We believe that the model could be useful in explaining biochemical energy recovery under cellular stress conditions.

Acknowledgments

The partial support of this work by a project of Vicerrectoría de Investigaciones of Universidad Nacional de Colombia (QUIPU No 201010013024) is gratefully acknowledged.

Supplementary Information

Supplementary information is available free of charge at http://jbcs.sbq.org.br as PDF file.

References

• ¹
  Winfree, A. T.; J. Chem. Educ. 1984, 61, 661.
• ²
  Tyson, J. In Frontiers in Mathematical Biology; Levin, S., ed.; Springer-Verlag Berlin Heidelberg: Berlin, Germany, 1994, p. 569-587.
• ³
  Pellitero, M. A.; Lamsfus, C. Á.; Borge, J.; J. Chem. Educ. 2012, 90, 82.
• ⁴
  Gray, P.; Griffiths J.; Combust. Flame 1989, 78, 87.
• ⁵
  Gray, P.; Ber. Bunsen-Ges. Phys. Chem. 1980, 84, 309.
• ⁶
  Griffiths, J. F.; Annu. Rev. Phys. Chem. 1985, 36, 77.
• ⁷
  Gray, P.; Kay, S. R.; Scott, S. K.; Proc. R. Soc. London, Ser. A 1988, 416, 321.
• ⁸
  Viel, F.; Jadot, F.; Bastin, G.; Automatica 1997, 33, 1437.
• ⁹
  Donner, J. S.; Thompson, S. A.; Kreuzer, M. P.; Baffou, G.; Quidant, R.; Nano Lett. 2012, 12, 2107.
• ¹⁰
  Gao, L.; Zhang, C.; Li, C.; Wang, L. V.; Appl. Phys. Lett. 2013, 102, 193705.
• ¹¹
  Sakaguchi, R.; Kiyonaka, S.; Mori, Y.; Curr. Opin. Biotechnol. 2015, 31, 57.
• ¹²
  Takei, Y.; Arai, S.; Murata, A.; Takabayashi, M.; Oyama, K.; Ishiwata, S.; Takeoka, S.; Suzuki, M.; ACS Nano 2013, 8, 198.
• ¹³
  Arai, S.; Lee, S.-C.; Zhai, D.; Suzuki, M.; Chang, Y. T.; Sci. Rep. 2014, 4, 6701.
• ¹⁴
  Tyson, J. In Computational Cell Biology; Fall, C.; Marland, E.; Wagner, J.; Tyson, J., eds.; Springer-Verlag New York: New York, USA, 2002, p. 230-260.
• ¹⁵
  Mankin, J. C.; Hudson, J. L.; Chem. Eng. Sci. 1984, 39, 1807.
• ¹⁶
  Luo, K.-M.; Lu, K.-T.; Hu, K.-H.; J. Loss Prev. Process Ind. 1997, 10, 141.
• ¹⁷
  Ball, R.; Ind. Eng. Chem. Res. 2013, 52, 922.
• ¹⁸
  Fujieda, S.; Kawahito, J.; Thermochim. Acta 1991, 183, 153.
• ¹⁹
  Lamprecht, I.; Schaarschmidt, B.; Plesser, T.; Thermochim. Acta 1987, 112, 95.
• ²⁰
  Nagy-Ungvárai, Z.; Kórös, E.; React. Kinet. Catal. Lett. 1985, 27, 83.
• ²¹
  Ágreda, J.; Barragán, D.; Gómez, A.; J. Therm. Anal. Calorim. 2003, 74, 875.
• ²²
  Vidal, C.; Noyau, A.; J. Am. Chem. Soc. 1980, 102, 6666.
• ²³
  Sahlin, K.; Sports Med. 2014, 44, 167.
• ²⁴
  MacIntosh, B. R.; Holash, R. J.; Renaud, J.-M.; J. Cell Sci. 2012, 125, 2105.
• ²⁵
  de Meis, L.; Bianconi, M. L.; Suzano, V. A.; FEBS Lett. 1997, 406, 201.
• ²⁶
  de Meis, L.; J. Biol. Chem. 2001, 276, 25078.
• ²⁷
  Nelson, M. I.; Sidhu, H. S.; J. Math. Chem. 2002, 31, 155.
• ²⁸
  The MathWorks, Inc.; MATLAB 2014b (version 8.4); Natick, Massachusetts, United States, 2014.
• ²⁹
  Gray, B. F.; Roberts, M. J.; Proc. R. Soc. London, Ser. A 1988, 416, 425.
• ³⁰
  Kay, S. R.; Scott, S. K.; Proc. R. Soc. London, Ser. A 1988, 416, 343.
• ³¹
  Gray, B. F.; Ball, R.; Proc. R. Soc. London, Ser. A 1999, 455, 163.
• ³²
  Zhang, J.; Zhou, L.; Ouyang, Q.; J. Phys. Chem. A 2007, 111, 1052.
• ³³
  Pullela, S. R.; Cristancho, D.; He, P.; Luo, D.; Hall, K. R.; Cheng, Z.; Phys. Chem. Chem. Phys. 2009, 11, 4236.
• ³⁴
  Allen, D. G.; Lamb, G. D.; Westerblad, H.; Physiol. Rev. 2008, 88, 287.
• ³⁵
  Ferguson, D. G.; Franzini-Armstrong, C.; Muscle Nerve 1988, 11, 561.
• ³⁶
  Higgins, J.; Proc. Natl. Acad. Sci. U. S. A. 1964, 51, 989.
• ³⁷
  Sel’Kov, E. E.; Eur. J. Biochem. 1968, 4, 79.
• ³⁸
  Roy, T.; Bhattacharjee, J. K.; Mallik, A. K.; Eur. Phys. J. E: Soft Matter Biol. Phys. 2011, 34, 19.
• ³⁹
  Mair, T.; Warnke, C.; Tsuji, K.; Müller, S. C.; Biophys. J. 2005, 88, 639.
• ⁴⁰
  Ørtenblad, N.; Macdonald, W. A.; Sahlin, K.; Biochem. J. 2009, 420, 161.
• ⁴¹
  Yang, J.-H.; Yang, L.; Qu, Z.; Weiss, J. N.; J. Biol. Chem. 2008, 283, 36321.
• ⁴²
  Kjelstrup, S.; Meis, L.; Bedeaux, D.; Simon, J.-M.; Eur. Biophys. J. 2008, 38, 59.
• ⁴³
  Bedeaux, D.; Kjelstrup, S.; Phys. Chem. Chem. Phys. 2008, 10, 7304.
• ⁴⁴
  Lervik, A.; Bresme, F.; Kjelstrup, S.; Bedeaux, D.; Miguel Rubi, J.; Phys. Chem. Chem. Phys. 2010, 12, 1610.
• ⁴⁵
  de Meis, L.; Biochem. Biophys. Res. Commun. 2000, 276, 35.
• ⁴⁶
  Ushimaru, M.; Fukushima, Y.; Biochem. Biophys. Res. Commun. 2007, 353, 799.
• ⁴⁷
  Muller, A. W. J.; Physiol. Chem. Phys. Med. NMR 1993, 25, 95.
• ⁴⁸
  Muller, A. W. J.; Biosystems 2005, 82, 93.
• ⁴⁹
  Muller, A. W. J.; Prog. Biophys. Mol. Biol. 1995, 63, 193.

Publication Dates

• Publication in this collection
  July 2018

History

• Received
  20 Sept 2017
• Accepted
  11 Jan 2018