Operational stability study of lactate biosensors: modeling, parameter identification, and stability analysis

Introduction: This paper investigates the operational stability of lactate biosensors, crucial devices in various biomedical and biotechnological applications. We detail the construction of an amperometric transducer tailored for lactate measurement and outline the experimental setup used for empirical validation.

Methods: The modeling framework incorporates Brown and Michaelis–Menten kinetics, integrating both distributed and discrete delays to capture the intricate dynamics of lactate sensing. To ascertain model parameters, we propose a nonlinear optimization method, leveraging initial approximations from the Brown model’s delay values for the subsequent model with discrete delays.

Results: Stability analysis forms a cornerstone of our investigation, centering on linearization around equilibrium states and scrutinizing the real parts of quasi-polynomials. Notably, our findings reveal that the discrete delay model manifests marginal stability, occupying a delicate balance between asymptotic stability and instability. We introduce criteria for verifying marginal stability based on characteristic quasi-polynomial roots, offering practical insights into system behavior.

Discussion: Qalitative examination of the model elucidates the influence of delay on dynamic behavior. We observe a transition from stable focus to limit cycle and period-doubling phenomena with increasing delay values, as evidenced by phase plots and bifurcation diagrams employing Poincaré sections. Additionally, we identify limitations in model applicability, notably the loss of solution positivity with growing delays, underscoring the necessity for cautious interpretation when employing delayed exponential function formulations. This comprehensive study provides valuable insights into the design and operational characteristics of lactate biosensors, offering a robust framework for understanding and optimizing their performance in diverse settings.

1 Introduction

1.1 Operational stability of biosensor vs. Lyapunov stability of the dynamic model

Operational stability of biosensors means “retention of activity of a protein or enzyme when in use” (Gibson, 1999).

It corresponds mainly with the same notion for dynamic systems from stability theory. Traditionally, enzyme–substrate interaction is simulated with the help of the Michaelis–Menten model, which is a nonlinear dynamic system. On the other hand, previous studies have paid little attention to the qualitative behavior of the model from the viewpoint of its stability. Partially, operational stability deals with the asymptotic nature of the stability notion, whereas biochemists conduct experiments during “finite time.”

In turn, the stability theory of dynamic systems, so-called Lyapunov stability, offers powerful tools for enhancing biochemical reactions modeling and saving the qualitative behavior of the systems. Moreover, each model is based on a series of assumptions about biochemical interactions, which allows us to check the validity of phenomena that cannot be verified experimentally, for example, mass action law with the distributed or discrete delays, which will be presented later.

1.2 Background for lactate measurement

Lactate characteristics: Lactate is an anion of lactic acid and is the final metabolite of the anaerobic breakdown of glucose. It is formed from pyruvate during the processes of glycolysis in the absence of oxygen (Figure 1) (Rabinowitz and Enerbäck, 2020) and is an important substance used in medicine as a marker of hypoxia and a number of other disorders, including diabetes and liver and kidney disorders (Pundir et al., 2016).

Figure 1

Figure 1. Scheme of reactions that underlie the functioning of a biosensor.

Lactate is also actively used in the food industry as an emulsifier, thickener, and acidity regulator. Various salts of lactic acid in the international classification of food additives are numbered E325, E326, E327, E328, and E329. Magnesium lactate is sometimes classified as an antioxidant (Standards, 2021). In addition, lactate can act as an indicator of the activity of bacteria during the fermentation process (Rawoof et al., 2020). It can indicate the freshness and quality of some products—wine (Gamella et al., 2010), juices (Trifirò et al., 1997), etc.

The concentration of lactate in products increases as they spoil due to changes in organoleptic properties. Due to the wide use of lactate as a marker of many processes in medicine and industry, accurate and effective methods of its diagnosis are necessary.

Lactate measurement methods: Traditionally, lactate is measured by colorimetry (Suman et al., 2005; Pundir et al., 2016), spectrophotometry (Pundir et al., 2016), fluorometry (Xue and Yeung, 1994), liquid chromatography (Biagi et al., 2012), and nuclear magnetic resonance (Hu et al., 2012; Park et al., 2015). These methods are certainly effective but have several disadvantages, such as high cost, length, and complexity of the analysis. Some of these methods also require complex preliminary sample preparation.

1.3 State-of-the-art lactate measurement with the help of biosensors

Biosensors offer effective and sensitive detection methods that can be used in medical institutions to measure the level of lactate in blood, sweat, and other biological fluids. In addition, new methods of lactate detection can be useful for enterprises that must control the processes of manufacturing food products and pharmaceuticals (Rattu et al., 2020).

Classification of lactate biosensors from the viewpoint of a bioselective element: A number of biosensor developments aimed at determining the concentration of lactate are known to be in development. Known biosensors can be divided according to the type of bioselective element: sensors based on lactate oxidase (LOx) and sensors based on lactate dehydrogenase (LHD) (Rassaei et al., 2013). In both cases, the substrate and product of the enzymatic reaction are lactate and pyruvate, respectively.

Basic reactions for lactate biosensors: There is a significant difference in these two reactions. For the LHD reaction, NAD+ is needed as a proton carrier in the reaction of dehydrogenation of lactate to pyruvate (Chatterjee et al., 2023). The lactate oxidase reaction can be much easier because, in it, the role of the proton acceptor is played by oxygen; that is, the use of this enzyme in biosensors does not require additional reagents (Lockridge et al., 1972).

Other variants of biosensors are based on a mixture of LHD and LOx. This configuration of the bioselective element, according to the results obtained by Chaubey et al. (2000), allows measuring lactate at lower concentrations than mono-enzyme biosensors, but this makes the analysis more expensive and is characterized by the complexity of manufacturing.

Classification of lactate biosensors from the viewpoint of the transducer: The main biosensor measurement methods used for lactate analysis are optical—electrochemiluminescence or fluorescence (Pundir and Narwal, 2017) and electrochemical—amperometric, potentiometric, or, less often, conductometric or impedimetric (Rathee et al., 2016). Moreover, various methods of improving the main characteristics of sensors have been applied to known biosensor systems—nanoparticles (Nesakumar et al., 2013; Azzouzi et al., 2015), other nanomaterials (Cui et al., 2007), complex, multi-stage methods of enzyme immobilization (Parra et al., 2006), and multi-enzyme membranes.

Although amperometric (Romero et al., 2010) and potentiometric (Mengarda et al., 2019) biosensors are usually monoenzymatic, that is, based on LHD or Lx, conductometric (Nguyen-Boisse et al., 2013) and impedimetric (Chan et al., 2017) biosensors usually use two-enzyme bioselective elements, for example, based on a mixture of lactate oxidase and peroxidase (Nguyen-Boisse et al., 2013) or a mixture of lactate dehydrogenase and pyruvate oxidase (Chan et al., 2017).

Motivation for lactate biosensor design: Therefore, in connection with the variety of works on the development of biosensor systems for measuring lactate and the prospect of its wide use in various areas of human life, we have concluded that the development of new biosensor methods for determining lactate is an urgent need.

Most existing biosensors are currently not ready for wide implementation and commercialization due to various limitations, such as insufficient sensitivity, selectivity, stability with respect to possible inferents, or too narrow a range of biosensor operation.

1.4 Michaelis–Menten model for enzyme kinetics

According to the model, an enzyme E combines with a substrate S to form an enzyme–substrate complex ES, characterized by a rate constant k₁. The resulting complex can dissociate into E and S (with a rate constant of k₋₁) or transform into a product P with a rate constant of k₂ (Berg et al., 2002).

The speed of the enzyme process is dependent on the ease of formation of the complex of the enzyme with the substrate. For low substrate concentrations, the reaction rate is proportional to the substrate concentration, while at higher concentrations, it tends toward a maximum value and becomes independent of substrate concentration. The general dependence of the rate of an enzyme reaction on substrate concentration is described by an equation called the Michaelis–Menten equation (Eq. 1):

$v=\frac{V_{\max }{n}_S}{K_M+n_S}.(1)$

Here, v is the rate of reaction (mol/s), V_max is the maximum reaction rate (mol/s), n_S is substrate concentration (mol/dm³), and K_M is the Michaelis–Menten constant (mol/dm³).

The Michaelis–Menten constant from Eq. 1 is an enzyme-specific quantity, dependent on substrate, temperature, and pH and independent of enzyme concentration. This constant is a measure of the affinity of the enzyme for the substrate. The lower the value of the constant, the higher the affinity of the enzyme for the substrate (Berg et al., 2002; Radomska, 2016).

1.5 Stability research on enzyme kinetics

Models of enzyme kinetics are based on compartmental systems, which are dynamic systems characterized by a network of interconnected nodes, each representing a reservoir or compartment where resources are stored (Blanchini et al., 2023). The system’s behavior is governed by the movement of resources, depicted as flows traveling along the edges connecting these compartments. Compartment-based dynamic systems serve as invaluable models across various disciplines, including physiologically based pharmacokinetics (Martsenyuk et al., 2012; Thompson and Beard, 2012), mathematical epidemiology (Reyné et al., 2022; Martsenyuk et al., 2021a; b), enzyme kinetics (Keener and Sneyd, 2009; Craciun et al., 2020; Martsenyuk et al., 2022), demography (Navarro Valencia et al., 2023), and ecology (Krishna et al., 2024). Stability analysis plays a pivotal role in understanding the behavior of these systems under different conditions (Blanchini et al., 2023). In recent years, researchers have made significant strides in advancing stability research methodologies, particularly in the context of compartmental systems with delay (Martsenyuk et al., 2013; Martsenyuk and Gandzyuk, 2013).

One common approach to stability analysis is linearization, which involves approximating nonlinear systems around equilibrium points. This technique has been extensively utilized in physiologically based pharmacokinetic models to assess the stability of drug distribution processes within the body (Martsenyuk et al., 2012). Lyapunov functions represent another powerful tool for stability analysis, offering a rigorous mathematical framework to prove the stability properties of compartment-based systems. In mathematical epidemiology, Lyapunov functions have been employed to establish global stability of disease-free and endemic equilibria in compartmental models of infectious diseases (Martsenyuk et al., 2021a).

1.6 Brief description of the work

Section 2 describes the materials, including the experiment and models used. The methods presented are related to parameter identification and stability research. Section 3 shows the results concerning the parameter identification for models using gamma-distributed delay and with two discrete delays. The qualitative analysis includes the existence and positiveness of the solutions, equilibrium states, marginal stability, and numerical research with the help of Poincaré sections. In Section 4, we discuss the results obtained and open problems.

The objective of the work is to offer the flowchart of the lactate biosensor design, including modeling, parameter identification, and stability analysis.

2 Materials and methods

2.1 Chemical compounds

The enzyme lactate oxidase obtained from Aerococcus viridans with an activity of 100 units (Sigma, United States) was used to create the biosensor. Bovine serum albumin (fraction V) (BSA) and a 25% aqueous solution of glutaraldehyde (GA), glycerol, and KCl were obtained from Sigma-Aldrich (United States). Stock 500 mM sodium L-lactate solution (Sigma-Aldrich, Switzerland) was used as a substrate. HEPES obtained from Sigma-Aldrich (United States) was used to prepare the buffer solution. Other inorganic compounds used in the work were of domestic production and had a purity level of “h.p.” and “p.d.a.”

2.2 Lactate biosensor design

The general scheme of amperometric transducers is shown in Figure 2. Platinum disc electrodes were used as amperometric transducers, manufactured in the laboratory of the Department of Biomolecular Electronics of the Institute of Biomolecular Biology and Geosciences using the following technology: a piece of platinum electrode with a diameter of 0.5 mm and a length of 3 mm was placed in a capillary tube with an outer diameter of 3.5 mm, and then the narrowed end of the capillary tube was sealed in a torch flame. The electrical connection between the platinum and the silver wire conductor was made by low-temperature soldering using Wood’s alloy. The open end of the capillary was filled with epoxy resin, with part of the conductor inside the capillary and part outside. A copper contact was soldered to the conductor to connect the transducer to the measuring unit.

Figure 2

Figure 2. General scheme of amperometric transducers. (A) Schematic view of the amperometric transducer. (B) Photograph of an amperometric transducer. (C) Scheme of the sensitive area of the transducer. (D) Photograph of the sensitive area of a transducer. 1–sensitive area, 2 –enzyme membrane, 3–platinum wire, 4–inner conductor (silver wire), electrical connection using 5–low-melting Wood’s alloy, 6–epoxy resin, and 7–contact area.

Amperometric measurements setup is shown in Figure 3. The PalmSens potentiostat (Palm Instruments BV, the Netherlands) was connected to an auxiliary platinum electrode, a silver chloride (Ag/AgCl) reference electrode, and working electrodes based on platinum disc electrodes. The potentiostat was connected to an 8-channel device (CH-8 multiplexer, Palm Instruments BV, the Netherlands), which allowed it to receive signals simultaneously from several working electrodes or biosensors (up to eight simultaneously). The distance between the auxiliary platinum electrode and all working biosensors during the measurement was the same (approximately 5 mm).

Figure 3

Figure 3. Scheme of the amperometric measurement setup. 1–auxiliary electrode, 2–measurement cell, 3–reference electrode, 4–working electrodes, 5–multiplexer, 6–potentiostat, and 7–PSTrace, measurement software for the potentiostat.

Preparation of bioselective membranes: Bioselective membranes were prepared by immobilizing proteins on the surface of a platinum disc electrode by covalent cross-linking of the enzyme in a bovine serum albumin matrix using glutaraldehyde as a cross-linking agent. The enzyme gel containing 5% lactotoxidase, 5% BSA, and 10% glycerol was mixed with a 1% glutaraldehyde solution in a 1:1 ratio, after which the mixture was applied to the sensitive area of the platinum disc electrode, and the membrane was air-dried for 20 min. After immobilization, the residual glutaraldehyde and unbound membrane components were washed off the membranes in a buffer solution for 5 min with constant stirring, changing the buffer several times.

2.3 Description of amperometric measurements

The measurement was carried out in an open measuring cell with a volume of 2 mL under constant stirring at room temperature. A fixed potential of +0.6 V was applied to the electrodes relative to the chloride silver reference electrode. The working buffer was 25 mM HEPES with pH 7.4. The required concentration of the substrate in the cell was set by adding aliquots of the standard solution, 500 mM sodium L-lactate, to the buffer.

The duration of one response (from the addition of the substrate to the signal output to the baseline) was approximately 4–5 min; between responses, the substrate was washed off the biosensor for 5 min, changing the buffer in the measuring cell several times. All measurements were performed in at least three replicates.

2.4 Models used

2.4.1 Model with continuously distributed delays

An application of delayed mass action law to enzyme kinetics was inspired by Brown’s model, formulated by Brown (1902), where complex C has a lifetime τ before being decayed. We called the reaction scheme

$\text{E}\left(\text{t}\right)+\text{S}\left(\text{t}\right)\stackrel{{\text{k}}_{\text{d}}}{\to }\text{E}\left(\text{t}+\tau \right)+\text{P}\left(\text{t}+\tau \right)$

an irreversible one-complex Brown’s (IR1CB) mechanism. In Martsenyuk et al. (2022), we offered the following model based on continuously distributed delays:

$\begin{array}{ll}\hfill \frac{d{n}_S\left(t\right)}{dt}& =-k_d{n}_E\left(t\right)n_S\left(t\right),\hfill \\ \hfill \frac{d{n}_E\left(t\right)}{dt}& =-k_d{n}_E\left(t\right)n_S\left(t\right)+k_d{\int }_{-{\tau }_M}^{0}f\left(s\right)n_E\left(t+s\right)n_S\left(t+s\right)ds,\hfill \\ \hfill \frac{d{n}_P\left(t\right)}{dt}& =k_d{\int }_{-{\tau }_M}^{0}f\left(s\right)n_E\left(t+s\right)n_S\left(t+s\right)ds,\hfill \end{array}(2)$

where for confidence level c ∈ (0, 1), we set ${\tau }_M≔{\tau }_{\mathit{min}}+\frac{m+1}{a}+\sqrt{\frac{(m+1)}{a^2(1-c)}}$, f(s) is the density function of the delay distribution, which was designed in the form of a gamma distribution:

$f\left(a,m,{\tau }_{\min },s\right)≔\left\{\begin{array}{cc}0\hfill & s\le {\tau }_{\min },\hfill \\ \frac{a^{m+1}}{\mathrm{\Gamma }\left(m+1\right)}{\left(s-{\tau }_{\min }\right)}^m{e}^{-a\left(s-{\tau }_{\min }\right)}\hfill & s>{\tau }_{\min },\hfill \end{array}\right.(3)$

where a, m, τ_min ≥ 0 are the parameters that determine the corresponding probability density function. Their roles and the ways of estimating were well-studied in Martsenyuk et al. (2022). The basic idea of Brown’s model shown in Model (2) does not include complex C directly but involves the model time τ required for complex forming–destroying. The model parameters were well-studied by Martsenyuk et al. (2022) and can be used as an initial approximation for the complex-based model in the next subsection.

2.4.2 Models with two discrete delays

The model extends the well-known irreversible one-complex Michaelis–Menten (IR1CMM) mechanism (Keener and Sneyd, 2009) (Section 1.4)

$\mathrm{E}+\mathrm{S}\underset{{\mathrm{k}}_{-1}}{\overset{{\mathrm{k}}_1}{\leftrightharpoons }}\mathrm{C}\stackrel{{\mathrm{k}}_2}{\to }\mathrm{E}+\mathrm{P}$

by adding the time durations τ₁ and τ₂ required for the entire fulfillment of the forward reactions. Mathematically, it corresponds to the time delays within dynamic systems. Hence, we consider the following set of elementary reactions:

$\begin{array}{ll}\hfill & \text{S}+\text{E}\stackrel{{\text{k}}_1,{\tau }_1}{\to }\text{C}\hfill \\ \hfill & \text{C}\stackrel{{\text{k}}_{-1}}{\to }\text{S}+\text{E}\hfill \\ \hfill & \text{C}\stackrel{{\text{k}}_2,{\tau }_2}{\to }\text{E}+\text{P},\hfill \end{array}$

which we call the irreversible one-complex with two delays Michaelis–Menthen (IR1C2DMM) mechanism. Based on the general approach described by Craciun et al. (2020), it yields the delayed model

$\begin{array}{ll}\hfill & \frac{d{n}_S}{dt}=k_{-1}{n}_C-k_1{n}_S\left(t-{\tau }_1\right)n_E\left(t-{\tau }_1\right)\hfill \\ \hfill & \frac{d{n}_E}{dt}=k_{-1}{n}_C+k_2{n}_C\left(t-{\tau }_2\right)-k_1{n}_S\left(t-{\tau }_1\right)n_E\left(t-{\tau }_1\right)\hfill \\ \hfill & \frac{d{n}_C}{dt}=k_1{n}_S\left(t-{\tau }_1\right)n_E\left(t-{\tau }_1\right)-k_2{n}_C\left(t-{\tau }_2\right)-k_{-1}{n}_C\hfill \\ \hfill & \frac{d{n}_P}{dt}=k_2{n}_C\left(t-{\tau }_2\right).\hfill \end{array}(4)$

For the solutions of Eq. 4, elements of which are the vector functions n_S, n_E, n_C, n_P ∈ C¹ ([−τ_max], 0], R⁴), we consider the following initial conditions:

$\begin{array}{ll}\hfill & n_S\left(t\right)={\widehat{n}}_S\left(t\right)\ge 0,n_E\left(t\right)={\widehat{n}}_E\left(t\right)\ge 0,n_C\left(t\right)={\widehat{n}}_C\left(t\right)\ge 0,n_P\left(t\right)={\widehat{n}}_P\left(t\right)\ge 0,t\in \left[-{\tau }_{\max },0\right),n_S\left(0\right)={\widehat{n}}_S^0>0,n_E\left(t\right)={\widehat{n}}_E^0>0,n_C\left(t\right)={\widehat{n}}_C^0>0,n_P\left(t\right)={\widehat{n}}_P^0>0.\hfill \end{array}(5)$

The value of n_P(t) can be found by direct integration, namely,

$n_P\left(t\right)=n_P\left(0\right)+k_2{\int }_0^t{n}_C\left(s-{\tau }_2\right)ds,t>0.(6)$

2.5 Methods

2.5.1 Parameter identification

As a result of the amperometric measurements described in Section 2, we obtain responses in the form of current I(t). We propose to use the relation of the current I(t) with n_P(t) as

$I\left(t\right)≔{\mathrm{\Lambda }}_m^0{n}_P\left(t\right)-K{\left(n_P\left(t\right)\right)}^{3/2},t>0.(7)$

In Martsenyuk et al. (2022), this relationship was evidenced for the specific conductance κ(t) of conductometric biosensors. Mathematical modeling of conductometric biosensors in terms of conductivity is presented in detail by Zouaoui et al. (2022). Provided fixed potential, we assume the linear dependence of the specific conductance and the current. So, we follow Eq. 7. Numerical modeling regarding amperometric biosensors is displayed by Simelevicius et al. (2012) and Hashem Zadeh et al. (2020).

In the following, we will denote product concentrations obtained as the solutions of the models as n_P,pred(t). In turn, the corresponding values of the current due to Eq. 7 will be I_pred(t). On the other hand, let I_exp(t) be the values of responses received as a result of the experiments.

The proposed parameter identification uses the currents I_exp,j (t_i) and I_pred,j (t_i), $j=\overline{1,m}$, experimentally and numerically obtained at the time instances t_i, $i=\overline{1,N}$¹ for given initial substrate concentrations n_S (0) = n_S,j.

Our goal is to estimate the parameters ${\mathrm{\Pi }}_{\text{IR}1\text{CB}}=\left\{k_d,a,m,{\tau }_{\mathit{min}},{\mathrm{\Lambda }}_m^0,K\right\}\in {\mathbf{R}}_{+}^6$, ${\mathrm{\Pi }}_{\text{IR}1\text{C}2\text{DMM}}=\left\{k_1,k_2,k_{-1},{\tau }_1,{\tau }_2,{\mathrm{\Lambda }}_m^0,K\right\}\in {\mathbf{R}}_{+}^7$ when applying Model (2) or (4), respectively.

We will estimate the parameters Π of the models solving the problems of optimization:

$\left.\begin{array}{cc}\hfill \mathtt{m}\mathtt{i}\mathtt{n}\mathtt{i}\mathtt{m}\mathtt{i}\mathtt{z}\mathtt{e}& \hfill J\left(\mathrm{\Pi }\right)\\ \mathtt{s}\mathtt{u}\mathtt{b}\mathtt{j}\mathtt{e}\mathtt{c}\mathtt{t}\mathtt{t}\mathtt{o}\hfill & g_i\left(\mathrm{\Pi }\right)\ge \mathrm{\Theta },i=\mathrm{1,2}.\hfill \end{array}\right\}(8)$

Here,

$J\left(\mathrm{\Pi }\right)≔{\left({\displaystyle \sum _{j=1}^m}\left({\displaystyle \sum _{i=1}^N}\left({\left(I_{\mathit{exp,j}}\left(t_i\right)-I_{\mathit{pred,j}}\left(t_i\right)\right)}^2\right)\right)\right)}^{1/2}(9)$

is the target function, and

$\begin{array}{cc}\hfill g_1\left(\mathrm{\Pi }\right)=\hfill & \mathrm{\Pi }-{\mathrm{\Pi }}_{\mathit{lower}}\ge \mathrm{\Theta },\hfill \\ \hfill g_2\left(\mathrm{\Pi }\right)=\hfill & {\mathrm{\Pi }}_{\mathit{upper}}-\mathrm{\Pi }\ge \mathrm{\Theta },\hfill \end{array}(10)$

are inequality constraints, where Θ is a null vector, and Π_lower and Π_upper are lower and upper bounds for the parameter values of corresponding dimensions.

An algorithm for the estimation of Π_IR1CB was described in detail by Martsenyuk et al. (2022). Here, we also apply it for the estimation of Π_IR1C2DMM. Moreover, the values of a, m, and τ_min, obtained as a result of the parameter identification in Eq, 2, will be used in order to set the initial estimate for τ₁ and τ₂ when estimating Eq. 4. Note that the mean value of time delay for Model (2) can be calculated as

$\mathtt{E}\left(\tau \right)≔{\tau }_{\mathit{min}}+\frac{m+1}{a}.$

So, we set the initial values of delays from Eq. 4 such that ${\tau }_1+{\tau }_2\approx \mathtt{E}(\tau )$.

2.5.2 Stability research using linearization technique

We will conduct research on local stability using the linearization technique. In this case, stability conditions are constructed based on a characteristic quasi-polynomial. The signs of the real parts of its roots are decisive for making conclusions about stability.

Namely, if all roots lie on the open left-half plane, then the equilibrium is locally asymptotically stable. If some of the roots have positive real parts, then the equilibrium is unstable. We focus our attention on the case when the roots lie on a closed left-half plane; that is, we have some simple roots lying on an imaginary axis. Such a situation is known as marginal stability. We distinguish this situation from the case of the pair of purely imaginary roots corresponding to periodic solutions.

When investigating the characteristic quasi-polynomial, a Padé approximant will be used in the form

$e^{-\lambda \tau }\approx \frac{1-\frac{\lambda \tau }{2}}{1+\frac{\lambda \tau }{2}},(11)$

allowing us to approximate the characteristic quasi-polynomial with the help of rational functions (Baker and Graves-Morris, 1996).

2.5.3 Bifurcation plots using Poincaré section

We use the Poincaré section technique to study the qualitative behavior of the models developed, which was primarily applied to the compartmental model by Martsenyuk et al. (2021a).

To begin, we gain a thorough understanding of Model (4), including its equations and parameters. We select parameters τ₁ and τ₂ to vary and the variables n_S, n_E, and n_C to focus on. Then, we simulate the system across different parameter values, generating time series data for the chosen variables.

Poincaré sections are constructed by intersecting the trajectory of the system with a defined plane n_E = d in the phase space, where d = (min_tn_E(t) + max_tn_E(t))/2. This section is determined by specific criteria; namely, we choose such points $(n_S^{\star },d,n_C^{\star })$ such that crossing a plane n_E = d will happen at sequential time instances t^⋆ + T, t^⋆ + 2T, t^⋆ + 3T, … , where T is a period value; that is, $n_S^{\star }=n_S(t^{\star })=n_S(t^{\star }+T)=n_S(t^{\star }+2T)=\dots$. We plot the sampled points in an (n_S, n_E)-space to visualize the Poincaré sections.

Repeating this process for various parameter values (τ₁, τ₂), we observe how the Poincaré sections change. Analyzing the patterns in the sections, we understand the system’s behavior, including the emergence of periodic orbits, chaotic dynamics, or transitions between different states.

Finally, we summarize the results by constructing a bifurcation plot combined with the corresponding n_S, n_E and (n_S, n_P) phase plots, showing how features of the Poincaré sections vary with parameter values. This comprehensive approach allows us to systematically explore the system’s behavior and gain insights into its dynamics.

3 Results

3.1 Parameter identification for a model using a gamma-distributed delay

The parameter identification technique mentioned in Section 2.5.1 was used for Model (2). The training data correspond to a set of time series of currents corresponding to given initial substrate concentrations n_S(0) equal to 0.1 mM, 0.5 mM, 1.0 mM, and 2.5 mM sequentially.

The goal is to estimate the parameters ${\mathrm{\Pi }}_{\text{IR}1\text{CB}}=\left\{k_d,a,m,{\tau }_{\mathit{min}},{\mathrm{\Lambda }}_m^0,K\right\}\in {\mathbf{R}}_{+}^6$. The initial values for the estimation were chosen as Π_init (see Table 1). The lower and upper bounds for Π_IR1CB are shown in Table 1 as Π_lower and Π_upper, respectively.

Table 1

Table 1. Parameter identification for Model (2).

The most valuable information obtained from Model (2) is the density of the distributed delay distribution (see Figure 4), which will be used in further estimations.

Figure 4

Figure 4. Density function given in 3 for the distributed delay τ value from Model Eq. (2).

The comparison of predicted and expected values of the currents for the optimal set of the parameters Π_opt is shown in Figure 5. We see that the parameter values, being optimal concerning the cost criterion (8), enable us to find the solution of Model (2) closest to the expected currents for the initial substrate concentration of 1.0 mM. For the smaller initial values of the substrate, we have predicted values smaller than the expected ones, whereas for those bigger than 1.0 mM, the predicted values are larger than the experimental ones. The explanation of such an effect lies in the special kind of stability of the model known as marginal stability, which will be evidenced further.

Figure 5

Figure 5. Comparison of the predicted and expected currents for Model (2) at the optimal values of the parameters Π_opt corresponding to given initial substrate concentrations n_S (0): (A) 0.1 mM, (B) 0.5 mM, (C) 1.0 mM, and (D) 2.5 mM.

3.2 Parameter identification for a model with two discrete delays

Model (4) requires the estimating the parameters ${\mathrm{\Pi }}_{\text{IR}1\text{C}2\text{DMM}}=\left\{k_1,k_2,k_{-1},{\tau }_1,{\tau }_2,{\mathrm{\Lambda }}_m^0,K\right\}\in {\mathbf{R}}_{+}^7$. They were obtained as a result of the solution of the optimization problem shown in Eqs 8–10. The training data described in Section 3.1 were used.

The initial values for the estimation were chosen as Π_init (see Table 2). The lower and upper bounds for Π_IR1C2DMM are shown in Table 2 as Π_lower and Π_upper, respectively.

Table 2

Table 2. Parameter identification for Model (4).

The comparison of predicted and expected values of the currents for the optimal set of the parameters Π_opt for Model (4) is shown in Figure 6. We see similar tendencies of expected and predicted values as we did for Model (2), with gamma-distributed delays.

Figure 6

Figure 6. Comparison of the predicted and expected currents for Model (4) at the optimal values of the parameters Π_opt corresponding to the given initial substrate concentrations n_S (0): (A) 0.1 mM, (B) 0.5 mM, (C) 1.0 mM, and (D) 2.5 mM.

3.3 Qualitative analysis

3.3.1 Existence and positiveness of the solutions

Given ${\widehat{n}}_S(t),{\widehat{n}}_E(t),{\widehat{n}}_C(t),{\widehat{n}}_P(t)\in {\mathbf{C}}^{+}[-{\tau }_{\max },0]$, because the right-hand sides of Eq. 4 imply the Lipschitz condition, there exists a unique trajectory of Eq. 4 starting from Eq. 5 (Hale and Lunel, 2013).

Henceforth, we will focus our attention on the positiveness of the solution of the system shown in Eq. 4.

The positiveness of n_P(t) follows directly from the positiveness of n_C(t). So, we prove the positiveness of n_S(t), n_E(t), and n_C(t) by contradiction.

3.3.1.1 Case without delays

First, we demonstrate the positiveness for Model (4) without delays, that is, if τ₁ = τ₂ = 0

$\begin{array}{ll}\hfill & \frac{d{n}_S}{dt}=k_{-1}{n}_C-k_1{n}_S{n}_E\hfill \\ \hfill & \frac{d{n}_E}{dt}=k_{-1}{n}_C+k_2{n}_C-k_1{n}_S{n}_E\hfill \\ \hfill & \frac{d{n}_C}{dt}=k_1{n}_S{n}_E-k_2{n}_C-k_{-1}{n}_C\hfill \\ \hfill & \frac{d{n}_P}{dt}=k_2{n}_C.\hfill \end{array}(12)$

Let us assume, for the sake of contradiction, that there is the smallest value among t_c, t_e, and t_s delivering non-positive solutions. Consider them sequentially.

Let t_c > 0 be the smallest instance of time such that n_C(t_c) = 0. From the first and second lines of Eq. 12, we get

$\begin{array}{ll}\hfill & \frac{d{n}_S}{dt}>-k_1{n}_S\left(t\right)n_E\left(t\right),\hfill \\ \hfill & \frac{d{n}_E}{dt}>-k_1{n}_S\left(t\right)n_E\left(t\right).\hfill \end{array}$

It implies that

$\begin{array}{ll}\hfill & n_S\left(t\right)>n_S\left(0\right)\exp \left(-k_1{\int }_0^t{n}_E\left(\mathit{\text{ξ}}\right)d\xi \right)>0,\hfill \\ \hfill & n_E\left(t\right)>n_E\left(0\right)\exp \left(-k_1{\int }_0^t{n}_E\left(\xi \right)d\xi \right)>0,\hfill \\ \hfill & t\in \left[0,t_c\right).\hfill \end{array}$

From the third line of Eq. 12, we have

$\frac{d{n}_C}{dt}>-k_2{n}_C\left(t\right)-k_{-1}{n}_C,t\in \left[0,t_c\right).$

Hence,

$n_C\left(t\right)>n_C\left(0\right)\exp \left(-k_2-k_{-1}\right)t,t\in \left[0,t_c\right).(13)$

By continuing Eq. 13, we have n_C(t_c) > 0, which contradicts the initial assumption.

Let t_e be the smallest instant that n_E(t_e) = 0. From the third part of Eq. 12, it follows that

$\frac{d{n}_C}{dt}>-k_2{n}_C\left(t\right)-k_{-1}{n}_C,$

and n_C(t) > 0 for t ∈ [0, t_e).

In turn, from the second part of Eq. 12, we have

$\frac{d{n}_E}{dt}>-k_1{n}_S\left(t\right)n_E\left(t\right),$

and

$n_E\left(t\right)>n_E\left(0\right)\exp \left(-k_1{\int }_0^t{n}_S\left(\xi \right)d\xi \right)>0$

for t ∈ [0, t_e],

Let t_s be the smallest instant that n_S (t_s) = 0. From the second line of Eq. 12, we get that n_C(t) > 0, t ∈ [0, t_s). Furthermore, from the first equation, we have

$\frac{d{n}_S}{dt}>-k_1{n}_S\left(t\right)n_E\left(t\right),$

and

$n_S\left(t\right)>n_S\left(0\right)\exp \left(-k_1{\int }_0^t{n}_E\left(\xi \right)d\xi \right)>0$

for t ∈ [0, t_s], which contradicts the assumption.

For the reasons given, we see that the solution of Eq. 12 exists and is positive for any positive initial values (n_S(0), n_E(0), n_C(0), n_P(0)) > 0.

3.3.1.2 Two discrete delays

It is natural to assume that the solutions of Eq. 4 are still positive for some sufficiently small τ₁ and τ₂. We aim to offer the conditions of positiveness.

The conditions will be based on the notion of the delayed exponential function. Given the value x ∈ R and delay τ > 0, the delayed exponential function is called (Angstmann et al., 2023)

$e_{\tau }^x≔{\displaystyle \sum _{n=0}^{+\infty }}\frac{{\left(x-n\tau \right)}^n}{\mathrm{\Gamma }\left(n+1\right)}\mathrm{\Theta }\left(\frac{x}{\tau }-n\right),$

where Θ(⋅) is the Heaviside function. This function has the property that $\frac{d}{dx}{e}_{\tau }^{\lambda x}=\lambda e_{\tau }^{\lambda (x-\tau )}$. Note that contrary to the “undelayed” exponential function, $e_{\tau }^x$ can also accept negative values.

To apply the approach of contradictions shown above, let t_c, t_e, and t_s be instances delivering non-positive solutions.

If t_c > 0 be the smallest instance that n_C (t_c) = 0, then from the first and second lines of Eq. 4

$\begin{array}{ll}\hfill & \frac{d{n}_S}{dt}>-k_1{n}_S\left(t-{\tau }_1\right)n_E\left(t-{\tau }_1\right),\hfill \\ \hfill & \frac{d{n}_E}{dt}>-k_1{n}_S\left(t-{\tau }_1\right)n_E\left(t-{\tau }_1\right).\hfill \end{array}$

Hence,

$\begin{array}{l}\hfill n_S\left(t\right)>n_S\left(0\right)e_{{\tau }_1}^{-k_1{n}_{E_0}t},\\ \hfill n_E\left(t\right)>n_{E_0}{e}_{{\tau }_1}^{-k_1{n}_S^{\max }t},\end{array}$

where $n_S^{\max }≔{\max }_u{n}_S(u)$.

From the third line of Eq. 4, we have

$\frac{d{n}_C}{dt}>-k_2{n}_C\left(t-{\tau }_2\right)-k_{-1}{n}_C,t\in \left[0,t_c\right).$

Consider t = t_c. Hence,

$\frac{d{n}_C\left(t_c\right)}{dt}>-k_2{n}_C\left(t_c-{\tau }_2\right),$

and we see that

$n_C\left(t\right)>n_C\left(0\right)e_{{\tau }_2}^{-k_2{n}_{E_0}t},t\in \left[0,t_c\right].$

Following the proof in the previous case, we conclude that the solutions of Eq. 4 are positive if rate parameters, initial conditions, and delays τ₁, τ₂ are such that for any t > 0

$e_{{\tau }_1}^{-k_1{n}_{E_0}t}>0,e_{{\tau }_1}^{-k_1{n}_S^{\max }t}>0,e_{{\tau }_2}^{-k_2{n}_{E_0}t}>0.$

3.3.2 Stability research

3.3.2.1 Equilibrium of the system

Let $({\bar{n}}_S,{\bar{n}}_E,{\bar{n}}_C,{\bar{n}}_P)$ be the equilibrium of Model (4). It should satisfy

$\begin{array}{ll}\hfill & k_{-1}{\bar{n}}_C-k_1{\bar{n}}_S{\bar{n}}_E=0,\hfill \\ \hfill & k_{-1}{\bar{n}}_C+k_2{\bar{n}}_C-k_1{\bar{n}}_S{\bar{n}}_E=0,\hfill \\ \hfill & k_1{\bar{n}}_S{\bar{n}}_E-k_2{\bar{n}}_C-k_{-1}{\bar{n}}_C=0,\hfill \\ \hfill & k_2{\bar{n}}_C=0.\hfill \end{array}$

It follows that

$\begin{array}{ll}\hfill & {\bar{n}}_C=0,\hfill \\ \hfill & k_1{\bar{n}}_S{\bar{n}}_E=0.\hfill \end{array}$

When analyzing Eq. 4, we see that there is a conserved quantity of enzyme because $\frac{d(n_E+n_C)}{dt}\equiv 0$. Hence, $n_E(t)+n_C(t)\equiv n_{E_0}$, where $n_{E_0}$ is the total amount of available enzyme. Thus, from the last equality of Section 3.3.2.1, we conclude that ${\bar{n}}_E=n_{E_0}$ and ${\bar{n}}_S=0$.

So, we have the unique equilibrium of Model (4), which is the substrate and complex free equilibrium (SCFE).

${\bar{n}}_S=0,{\bar{n}}_E=n_{E_0},{\bar{n}}_C=0,{\bar{n}}_P=\text{const}.$

The value of ${\bar{n}}_P$ is bounded and can be determined by Eq. 6. Note that it is related to the initial values of both n_S and n_E. n_P is not included in the right-hand sides of Eq. 4.

3.3.2.2 Marginal stability

The Jacobian matrix at SCFE for Model (4) is given by

$\mathbf{J}={\left[\begin{array}{cccc}\hfill -k_1{\bar{n}}_E{e}^{-\lambda {\tau }_1}\hfill & \hfill -k_1{\bar{n}}_S{e}^{-\lambda {\tau }_1}\hfill & \hfill k_{-1}\hfill & \hfill 0\hfill \\ \hfill -k_1{\bar{n}}_E{e}^{-\lambda {\tau }_1}\hfill & \hfill -k_1{\bar{n}}_S{e}^{-\lambda {\tau }_1}\hfill & \hfill k_{-1}+k_2{e}^{-\lambda {\tau }_2}\hfill & \hfill 0\hfill \\ \hfill k_1{\bar{n}}_E{e}^{-\lambda {\tau }_1}\hfill & \hfill k_1{\bar{n}}_S{e}^{-\lambda {\tau }_1}\hfill & \hfill -k_2{e}^{-\lambda {\tau }_2}-k_1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill k_2{e}^{-\lambda {\tau }_2}\hfill & \hfill 0\hfill \end{array}\right]}_{\text{SCFE}}=\left[\begin{array}{cccc}\hfill -k_1{n}_{E_0}{e}^{-\lambda {\tau }_1}\hfill & \hfill 0\hfill & \hfill k_{-1}\hfill & \hfill 0\hfill \\ \hfill -k_1{n}_{E_0}{e}^{-\lambda {\tau }_1}\hfill & \hfill 0\hfill & \hfill k_{-1}+k_2{e}^{-\lambda {\tau }_2}\hfill & \hfill 0\hfill \\ \hfill k_1{n}_{E_0}{e}^{-\lambda {\tau }_1}\hfill & \hfill 0\hfill & \hfill -k_2{e}^{-\lambda {\tau }_2}-k_1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill k_2{e}^{-\lambda {\tau }_2}\hfill & \hfill 0\hfill \end{array}\right].$

The characteristic quasi-polynomial for $({\bar{n}}_S,{\bar{n}}_E,{\bar{n}}_C,{\bar{n}}_P)$ is given by

$\chi \left(\lambda \right)={\lambda }^4+k_{-1}{\lambda }^3+k_1\left({\bar{n}}_S+{\bar{n}}_E\right){\lambda }^3{e}^{-\lambda {\tau }_1}+k_2{\lambda }^3{e}^{-\lambda {\tau }_2}+k_1{k}_2{\bar{n}}_E{\lambda }^2{e}^{-\lambda \left({\tau }_1+{\tau }_2\right)}={\lambda }^2{\chi }_1\left(\lambda \right),$

where

${\chi }_1\left(\lambda \right)={\lambda }^2+k_{-1}\lambda +k_1\left({\bar{n}}_S+{\bar{n}}_E\right)\lambda e^{-\lambda {\tau }_1}+k_2\lambda e^{-\lambda {\tau }_2}+k_1{k}_2{\bar{n}}_E{e}^{-\lambda \left({\tau }_1+{\tau }_2\right)}.$

Hence, stability analysis can be reduced to obtaining conditions of non-positive values of real parts of the quasi-polynomial of χ₁(λ).

For SCFE, we have the following quasi-polynomial:

${\chi }_1\left(\lambda \right)={\lambda }^2+k_{-1}\lambda +k_1{n}_{E_0}\lambda e^{-\lambda {\tau }_1}+k_2\lambda e^{-\lambda {\tau }_2}+k_1{k}_2{\bar{n}}_E{e}^{-\lambda \left({\tau }_1+{\tau }_2\right)}.$

When applying the Padé approximation shown in Eq. 11 to the characteristic quasi-polynomial, we get

$\begin{array}{ll}\hfill {\chi }_1\left(\lambda \right)\approx & \left(k_1{\bar{n}}_E{\lambda }^3{\tau }_1^2{k}_2{\tau }_2-2k_1{\bar{n}}_E{\lambda }^3\tau 1^2-k_1{\bar{n}}_E{\lambda }^4{\tau }_1^2{\tau }_2+-2k_1{\bar{n}}_E{\lambda }^2\tau 1^2{k}_2+4k_1{\bar{n}}_E{\lambda }^2{\tau }_2\right.\hfill \\ \hfill & -4k_1{\bar{n}}_E\lambda k_2{\tau }_2+8k_1{\bar{n}}_E\lambda +8k_1{\bar{n}}_E{k}_2+{\lambda }^5{\tau }_1^2{\tau }_2+{\lambda }^4\tau 1^2{k}_{-1}{\tau }_2+2{\lambda }^4{\tau }_1^2\hfill \\ \hfill & -{\lambda }^4{\tau }_1^2{k}_2{\tau }_2+4{\lambda }^4{\tau }_1{\tau }_2+2{\lambda }^3{\tau }_1^2{k}_{-1}+2{\lambda }^3{\tau }_1^2{k}_2+4{\lambda }^3{\tau }_1{k}_{-1}{\tau }_2-4{\lambda }^3{\tau }_1{k}_2{\tau }_2\hfill \\ \hfill & \left.+8{\lambda }^3{\tau }_1+4{\lambda }^3{\tau }_2+8{\lambda }^2{\tau }_1{k}_{-1}+8{\lambda }^2{\tau }_1{k}_2+4{\lambda }^2{k}_{-1}{\tau }_2-4{\lambda }^2{k}_2{\tau }_2+8{\lambda }^2+8\lambda k_{-1}+8\lambda k_2\right)/\left({\lambda }^3{\tau }_1^2{\tau }_2+2{\lambda }^2{\tau }_1^2+4{\lambda }^2{\tau }_1{\tau }_2+8\lambda {\tau }_1+4\lambda {\tau }_2+8\right).\hfill \end{array}$

Setting the values for ${\bar{n}}_E=1.139071mM$, k₁, k₂, k₋₁ from Π_opt (Table 2), we obtain

${\chi }_1\left(\lambda \right)=0.0992985372e^{-\lambda {\tau }_1}\lambda +0.0121224945e^{-\lambda {\tau }_1}{e}^{-\lambda {\tau }_2}+{\lambda }^2+0.1220813\lambda e^{-\lambda {\tau }_2}+0.1048326\lambda .$

In turn, applying the Padé approximation yields

$\begin{array}{ll}\hfill {\chi }_1\left(\lambda \right)\approx & \left(8\left(\left({\lambda }^5{\tau }_1^2{\tau }_2\right)/8+\left(-0.9323778976{\lambda }^4{\tau }_1^2{\tau }_2\right)/64+\left({\lambda }^4{\tau }_1^2\right)/4+\left({\lambda }^4{\tau }_1{\tau }_2\right)/2\right.\right.+\left(0.1939599121{\lambda }^3{\tau }_1^2{\tau }_2\right)/128+\left(0.2552307255{\lambda }^3{\tau }_1^2\right)/8+\left(-4.4156672{\lambda }^3{\tau }_1{\tau }_2\right)/512\hfill \\ \hfill & +{\lambda }^3{\tau }_1+\left({\lambda }^3{\tau }_2\right)/2+\left(-0.048489978{\lambda }^2{\tau }_1^2\right)/16+\left(1.8153112{\lambda }^2{\tau }_1\right)/8+\left(0.3281993488{\lambda }^2{\tau }_2\right)/8\left.\left.+{\lambda }^2+\left(-0.048489978\lambda {\tau }_2\right)/8+0.3262124372\lambda +0.0121224945\right)\right)/\left({\lambda }^3{\tau }_1^2{\tau }_2+2{\lambda }^2{\tau }_1^2\right.\hfill \\ \hfill & \left.+4{\lambda }^2{\tau }_1{\tau }_2+8\lambda {\tau }_1+4\lambda {\tau }_2+8\right).\hfill \end{array}$

Consider two special cases of Eq. 4:

Case 1 without delays, that is, τ₁ = τ₂ = 0. Then,

${\chi }_1\left(\lambda \right)={\lambda }^2+0.3262124372\lambda +0.0121224945,$

and we have the roots of χ₁(λ):

${\lambda }_1=-0.28344384,{\lambda }_2=-0.04276859,$

which means marginal stability in this case.

Case 2 with optimal values of the delays from Π_opt (Table 2), that is, τ₁ = 15.69849 and τ₂ = 5.357829. It holds

$\begin{array}{ll}\hfill {\chi }_1\left(\lambda \right)=& 0.0992985372e^{-15.69849\lambda }\lambda +0.0121224945e^{-15.69849\lambda }{e}^{-5.357829\lambda }\hfill \\ \hfill & +{\lambda }^2+0.1220813\lambda e^{-5.357829\lambda }+0.1048326\lambda .\hfill \end{array}$

3.3.3 Effect of delays on qualitative behavior

We investigated the qualitative behavior of the model by selecting specific parameter values, denoted as Π_opt, and examining the influence of the delay τ₁ on Model 4’s dynamics, provided that τ₂ = 5.357829 is fixed. To explore this influence, we varied the parameter τ₁ over the interval [0, 35].

Figure 7 illustrates the construction of Poincaré sections in 3D plots for various values of τ₁. These sections provide insights into the system’s behavior at specific points in its phase space. As it can be seen for the values of τ₁ and τ₂ greater than in Π_opt, we lose the positivity of the solutions. This phenomenon was evidenced in Section 3.3.1.2, where it was shown that in contrast to the model without delay, positivity requires restricting the parameters in Eq. (4).

Figure 7

Figure 7. Constructing Poincaré sections based on trajectories (n_S(t), n_E(t), n_C(t)) of Model (4) with the parameters Π_opt in dependence of (A) τ₁ =15.69849, (B) τ₁ = 25, (C) τ₁ = 32, and (D) τ₁ = 35. Poincaré sections (in red) are obtained as a result of crossing with the plane n_E = d, where d =(min_tn_E(t) + max_tn_E(t))/2.

Furthermore, Figure 8 presents the corresponding Poincaré sections plotted in the (n_S, n_E)-plane, offering a clearer visualization of the system’s trajectory crossings. A thorough inspection of the sections says that for the initial trajectory (n_S(t), n_E(t), n_C(t)) behaves at the equilibrium state as a stable node. Then, at less than τ ≈ 25, Hopf bifurcation appears, and the limit cycle starts from a small radius of approximately 10^–4 and is extended sequentially to 10^–2. Simultaneously, period doubling occurs.

Figure 8

Figure 8. Poincaré sections (in blue) in the (n_S, n_E)-plane of the model (4) with the parameters Π_opt in dependence of (A) τ₁ = 15.69849, (B) τ₁ = 25, (C) τ₁ = 32, and (D) τ₁ = 35.

In Figure 9, we depict the bifurcation diagram for the variable n_S resulting from the parameter variation of τ₁. It reveals that the equilibrium state behaves as a stable node within the range of τ₁ values from 0 to 25. Subsequently, a bifurcation occurs, leading to a transition to a limit cycle. This period-doubling phenomenon is further supported by Figure 10, which displays the (n_S, n_E)-phase plane for different initial states.

Figure 9

Figure 9. Bifurcation diagram of n_S with respect to the parameter τ₁.

Figure 10

Figure 10. (n_S, n_E)-phase plots of Model (4) with the parameters Π_opt in dependence of (A) τ₁=15.69849, (B) τ₁=25, (C) τ₁=32, and (D) τ₁=35. Trajectories are constructed at different initial values of n_S (0):–0.1 mM, –0.5 mM, –1.0 mM, and –2.5 mM. Starting points are indicated in blue.

For τ₁ values exceeding 32, the system’s behavior appears likely to become chaotic. This transition is occurring through period doubling, although further investigation is warranted. It is noteworthy that the system exhibits marginal stability, as evident from Figure 11. This figure illustrates how different initial values of n_S (0) lead the solution n_P to converge toward either a stable node or a limited cycle.

Figure 11

Figure 11. (n_S, n_P)-phase plots of Model (4) with the parameters Π_opt in dependence of (A) τ₁=15.69849, (B) τ₁=25, (C) τ₁=32, and (D) τ₁=35. Trajectories are constructed at different initial values of n_S (0):–0.1 mM, –0.5 mM, –1.0 mM, and –2.5 mM. Starting points are indicated in blue.

Special attention should be paid to the positivity of the solutions, as we observed from the numerical modeling.

4 Discussion

The use of gamma-distributed delays and discretely distributed delay models is simultaneously important because the parameter values from the IR1CB model can be used for the IR1C2DMM model. The gamma-distributed delays model allows us to estimate the mean value of the delay easily.

Moreover, based on the results of numerical experiments, we can conclude that for the optimal values of the parameters Π_opt (Tables 1, 2), we have that τ₁ + τ₂ ≈ E(τ), where τ₁, τ₂ are the discrete delays from Eq. 4, and τ is the distributed delay from Eq. 2 (see also Figure 4 for delay density distribution).

The use of the Levenberg–Marquardt algorithm is essentially determined by choosing the initial approximation of the parameter values. It can be improved by applying AI techniques and constructing and tuning neural networks of the appropriate architecture for future consideration.

Artificial intelligence (AI) models like recurrent NNs can also be used to model enzyme–substrate interaction (we will try this in the future). The advantage of an AI model is that it makes more accurate predictions. The disadvantages of using AI are (1) the black box problem, meaning that the models are not based on any biochemical assumptions but only neural network architecture and (2) the overfitting problem, meaning that when we try to fit the outputs as accurately as we can, we may not see the general tendencies in reactions.

The system is characterized by marginal stability, which is between Lyapunov stability and instability. It corresponds to the objective of an electrochemical biosensor as a measuring device. Each initial state, including substrate concentration changes, has its “own” concentration of the product at an infinite time. We can reformulate it as the definition of operating stability.

The question arises of how to apply and interpret the marginal stability condition. As was shown, the marginal stability condition can be reduced by checking the real parts of the polynomial roots. We see that the roots depend on the reaction rate parameters (k₁, k₂, k₋₁), time delays τ₁, τ₂, and the initial concentration of $n_{E_0}$.

5 Conclusion

The equilibrium of the system was demonstrated in this article. From the viewpoint of enzymatic reactions and based on the proposed model, we saw that the system would be in equilibrium in the case of the absence of substrate and complex, that is, ${\bar{n}}_S\equiv 0$, n_C ≡ 0, as well as the concentration of the enzyme equal to the total amount of available enzyme. The concentration of the product should be a constant value calculated in accordance with Eq. 6 as the area under curve n_C(t), that is, ${\bar{n}}_P\equiv const$. In addition, ${\bar{n}}_P$ is determined by the initial value n_S (0).

In the work, we also considered a model with two discrete delays where τ₁ and τ₂ are time durations required for the entire fulfillment of forwarded reactions. To be precise, τ₁ is the time needed for the substrate to bind with the enzyme (time needed for complex formation), and τ₂ is the time needed for the complex to break down into enzyme and product. It was numerically shown that increasing those two parameters could result in a loss of stability. Suppose we have optimal values ${\tau }_1^0$ and ${\tau }_2^0$ and they correspond to a stable solution. If ${\tau }_1>{\tau }_1^0$, more time is required for forming the complex at $\text{S}+\text{E}\stackrel{{\text{k}}_1\text{,}{\tau }_1}{\to }\text{C}$. In other words, less complex will be produced compared with the optimal value. In turn, if ${\tau }_2>{\tau }_2^0$, this means that more time will be needed for the $\text{C}\stackrel{{\text{k}}_1\text{,}{\tau }_2}{\to }\text{E}+\text{P}$ reaction, leading to less enzyme and product than in the steady state. We remember that in the steady state, it has been proven that we have no complex, and the enzyme content is maximal. On the other hand, an increase in τ₂ results in having more complex but less product. This will affect the final content of the product, that is, the biosensor index.

It can also be noted that no unique parameters are set for fitting the model to the expected data at different values (Figures 5, 6), as the rates are not constant but can be functions of the substance concentrations. In particular, the times required for the reactions $\text{S}+\text{E}\stackrel{{\text{k}}_1\text{,}{\tau }_1}{\to }\text{C}$ and $\text{C}\stackrel{{\text{k}}_1\text{,}{\tau }_2}{\to }\text{E}+\text{P}$ are not constant. They likely depend on the concentrations of S, E (for τ₁), and C (for τ₂). These data should be determined experimentally. Moreover, the models can be extended by involving other complexes created by reactions with other biosensor auxiliary components. For example, BSA is known to create a complex with an enzyme known as CLEA (Shah et al., 2006). In other words, incorporating additional variables and increasing the system dimension can improve the model fitting.

These findings provide valuable insights into the system’s dynamics and highlight the importance of parameter exploration in understanding the qualitative behavior of enzyme kinetics models.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

VM: writing–review and editing, writing–original draft, visualization, validation, supervision, software, methodology, investigation, and conceptualization. OS: writing–review and editing, writing–original draft, validation, supervision, resources, project administration, methodology, funding acquisition, formal analysis, data curation, and conceptualization. AK-W: writing–review and editing, writing–original draft, visualization, validation, methodology, investigation, and formal analysis. AS: writing–review and editing, writing–original draft, validation, project administration, methodology, investigation, funding acquisition, data curation, and conceptualization. KB: writing–review and editing, writing–original draft, resources, methodology, investigation, and data curation.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The work was supported by the National Research Foundation of Ukraine in the framework of the competition of projects for research and development, “Science for the reconstruction of Ukraine in the war and post-war periods” (project 2022.01/0043).

Acknowledgments

The authors are grateful to the reviewers for valuable comments that allowed the paper to be significantly improved. All activities related to preparing and publishing the paper were financed by the University of Bielsko-Biala within the framework of the discipline “Technical Informatics and Telecommunication.”

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Without loss of generality, we assume that all time series have the same time span.

References

Angstmann, C. N., Burney, S.-J. M., Henry, B. I., Jacobs, B. A., and Xu, Z. (2023). A systematic approach to delay functions. Mathematics 11, 4526. doi:10.3390/math11214526

CrossRef Full Text | Google Scholar

Azzouzi, S., Rotariu, L., Benito, A. M., Maser, W. K., Ali, M. B., and Bala, C. (2015). A novel amperometric biosensor based on gold nanoparticles anchored on reduced graphene oxide for sensitive detection of l-lactate tumor biomarker. Biosens. Bioelectron. 69, 280–286. doi:10.1016/j.bios.2015.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Baker, G. A., and Graves-Morris, P. (1996). Padé approximants. Second Edition. Cambridge: Cambridge University Press. doi:10.1017/cbo9780511530074

CrossRef Full Text | Google Scholar

Berg, J., Tymoczko, J., and Stryer, L. (2002). Biochemistry. Fifth Edition. Germany: W. H. Freeman.

Google Scholar

Biagi, S., Ghimenti, S., Onor, M., and Bramanti, E. (2012). Simultaneous determination of lactate and pyruvate in human sweat using reversed-phase high-performance liquid chromatography: a noninvasive approach. Biomed. Chromatogr. 26, 1408–1415. doi:10.1002/bmc.2713

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanchini, F., Bolzern, P., Colaneri, P., De Nicolao, G., and Giordano, G. (2023). Optimal control of compartmental models: the exact solution. Automatica 147, 110680. doi:10.1016/j.automatica.2022.110680

CrossRef Full Text | Google Scholar

Brown, A. J. (1902). Xxxvi.—enzyme action. J. Chem. Soc. Trans. 81, 373–388. doi:10.1039/ct9028100373

CrossRef Full Text | Google Scholar

Chan, D., Barsan, M. M., Korpan, Y., and Brett, C. M. (2017). L-lactate selective impedimetric bienzymatic biosensor based on lactate dehydrogenase and pyruvate oxidase. Electrochimica Acta 231, 209–215. doi:10.1016/j.electacta.2017.02.050

CrossRef Full Text | Google Scholar

Chatterjee, Y., Bhowal, B., Gupta, K. J., Pareek, A., and Singla-Pareek, S. L. (2023). Lactate dehydrogenase superfamily in rice and arabidopsis: understanding the molecular evolution and structural diversity. Int. J. Mol. Sci. 24, 5900. doi:10.3390/ijms24065900

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaubey, A., Pande, K., Singh, V., and Malhotra, B. (2000). Co-immobilization of lactate oxidase and lactate dehydrogenase on conducting polyaniline films. Anal. Chim. Acta 407, 97–103. doi:10.1016/s0003-2670(99)00797-7

CrossRef Full Text | Google Scholar

Craciun, G., Mincheva, M., Pantea, C., and Yu, P. Y. (2020). Delay stability of reaction systems. Math. Biosci. 326, 108387. doi:10.1016/j.mbs.2020.108387

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, X., Li, C. M., Zang, J., and Yu, S. (2007). Highly sensitive lactate biosensor by engineering chitosan/PVI-os/CNT/LOD network nanocomposite. Biosens. Bioelectron. 22, 3288–3292. doi:10.1016/j.bios.2007.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Gamella, M., Campuzano, S., Conzuelo, F., Curiel, J., Muñoz, R., Reviejo, A., et al. (2010). Integrated multienzyme electrochemical biosensors for monitoring malolactic fermentation in wines. Talanta 81, 925–933. doi:10.1016/j.talanta.2010.01.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibson, T. D. (1999). Biosensors: the stabilité problem. Analusis 27, 630–638. doi:10.1051/analusis:1999270630

CrossRef Full Text | Google Scholar

Hale, J. K., and Lunel, S. M. V. (2013). Introduction to functional differential equations, 99. Germany: Springer Science and Business Media.

Google Scholar

Hashem Zadeh, S. M., Heidarshenas, M., Ghalambaz, M., Noghrehabadi, A., and Saffari Pour, M. (2020). Numerical modeling and investigation of amperometric biosensors with perforated membranes. Sensors 20, 2910. doi:10.3390/s20102910

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, S., Yoshihara, H. A., Bok, R., Zhou, J., Zhu, M., Kurhanewicz, J., et al. (2012). Use of hyperpolarized [1-13c]pyruvate and [2-13c]pyruvate to probe the effects of the anticancer agent dichloroacetate on mitochondrial metabolism in vivo in the normal rat. Magn. Reson. Imaging 30, 1367–1372. doi:10.1016/j.mri.2012.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

J. Keener, and J. Sneyd (2009). Mathematical physiology (New York: Springer). doi:10.1007/978-0-387-75847-3

CrossRef Full Text | Google Scholar

Krishna, S., Peterson, V., Listmann, L., and Hinners, J. (2024). Interactive effects of viral lysis and warming in a coastal ocean identified from an idealized ecosystem model. Ecol. Model. 487, 110550. doi:10.1016/j.ecolmodel.2023.110550

CrossRef Full Text | Google Scholar

Lockridge, O., Massey, V., and Sullivan, P. A. (1972). Mechanism of action of the flavoenzyme lactate oxidase. J. Biol. Chem. 247, 8097–8106. doi:10.1016/s0021-9258(20)81814-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Martsenyuk, V., Augustynek, K., and Urbas, A. (2021a). On qualitative analysis of the nonstationary delayed model of coexistence of two-strain virus: stability, bifurcation, and transition to chaos. Int. J. Non-Linear Mech. 128, 103630. doi:10.1016/j.ijnonlinmec.2020.103630

PubMed Abstract | CrossRef Full Text | Google Scholar

Martsenyuk, V., Bernas, M., and Klos-Witkowska, A. (2021b). Two-strain covid-19 model using delayed dynamic system and big data. IEEE Access 9, 113866–113878. doi:10.1109/access.2021.3104519

CrossRef Full Text | Google Scholar

Martsenyuk, V., Klos-Witkowska, A., Dzyadevych, S., and Sverstiuk, A. (2022). Nonlinear analytics for electrochemical biosensor design using enzyme aggregates and delayed mass action. Sensors 22, 980. doi:10.3390/s22030980

PubMed Abstract | CrossRef Full Text | Google Scholar

Martsenyuk, V. P., Andrushchak, I. E., and Melenchuk, I. B. (2012). Method of construction and determination of approximate solutions of the model of pharmacokinetics of nanoparticles. J. Automation Inf. Sci. 44, 32–43. doi:10.1615/jautomatinfscien.v44.i8.40

CrossRef Full Text | Google Scholar

Martsenyuk, V. P., Andrushchak, I. Y., and Gandzyuk, N. M. (2013). Constructing exponential estimates in compartmental systems with distributed delays: an approach based on the hale–lunel inequality. Cybern. Syst. Analysis 49, 347–352. doi:10.1007/s10559-013-9517-0

CrossRef Full Text | Google Scholar

Martsenyuk, V. P., and Gandzyuk, N. M. (2013). Stability estimation method for compartmental models with delay. Cybern. Syst. Analysis 49, 81–85. doi:10.1007/s10559-013-9488-1

CrossRef Full Text | Google Scholar

Mengarda, P., Dias, F. A., Peixoto, J. V., Osiecki, R., Bergamini, M. F., and Marcolino-Junior, L. H. (2019). Determination of lactate levels in biological fluids using a disposable ion-selective potentiometric sensor based on polypyrrole films. Sensors Actuators B Chem. 296, 126663. doi:10.1016/j.snb.2019.126663

CrossRef Full Text | Google Scholar

Navarro Valencia, V. A., Díaz, Y., Pascale, J. M., Boni, M. F., and Sanchez-Galan, J. E. (2023). Using compartmental models and particle swarm optimization to assess dengue basic reproduction number r0 for the republic of Panama in the 1999-2022 period. Heliyon 9, e15424. doi:10.1016/j.heliyon.2023.e15424

PubMed Abstract | CrossRef Full Text | Google Scholar

Nesakumar, N., Sethuraman, S., Krishnan, U. M., and Rayappan, J. B. B. (2013). Fabrication of lactate biosensor based on lactate dehydrogenase immobilized on cerium oxide nanoparticles. J. Colloid Interface Sci. 410, 158–164. doi:10.1016/j.jcis.2013.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen-Boisse, T. T., Saulnier, J., Jaffrezic-Renault, N., and Lagarde, F. (2013). Highly sensitive conductometric biosensors for total lactate, d- and l-lactate determination in dairy products. Sensors Actuators B Chem. 179, 232–239. doi:10.1016/j.snb.2012.10.021

CrossRef Full Text | Google Scholar

Park, J. M., Josan, S., Mayer, D., Hurd, R. E., Chung, Y., Bendahan, D., et al. (2015). Hyperpolarized 13c NMR observation of lactate kinetics in skeletal muscle. J. Exp. Biol. 218, 3308–3318. doi:10.1242/jeb.123141

PubMed Abstract | CrossRef Full Text | Google Scholar

Parra, A., Casero, E., Vázquez, L., Pariente, F., and Lorenzo, E. (2006). Design and characterization of a lactate biosensor based on immobilized lactate oxidase onto gold surfaces. Anal. Chim. Acta 555, 308–315. doi:10.1016/j.aca.2005.09.025

CrossRef Full Text | Google Scholar

Pundir, C. S., and Narwal, V. (2017). A bird’s eye view of lactate biosensors. J. Intensive Crit. Care 03. doi:10.21767/2471-8505.100089

CrossRef Full Text | Google Scholar

Pundir, C. S., Narwal, V., and Batra, B. (2016). Determination of lactic acid with special emphasis on biosensing methods: a review. Biosens. Bioelectron. 86, 777–790. doi:10.1016/j.bios.2016.07.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Rabinowitz, J. D., and Enerbäck, S. (2020). Lactate: the ugly duckling of energy metabolism. Nat. Metab. 2, 566–571. doi:10.1038/s42255-020-0243-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Radomska (2016). Praca dokt_mradomska.pdf. Available at: https://depotuw.ceon.pl/bitstream/handle/item/2777/praca%20dokt_%20MRadomska.pdf?sequence=1 (Accessed on April 11, 2023).

Google Scholar

Rassaei, L., Olthuis, W., Tsujimura, S., Sudhölter, E. J. R., and van den Berg, A. (2013). Lactate biosensors: current status and outlook. Anal. Bioanal. Chem. 406, 123–137. doi:10.1007/s00216-013-7307-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rathee, K., Dhull, V., Dhull, R., and Singh, S. (2016). Biosensors based on electrochemical lactate detection: a comprehensive review. Biochem. Biophysics Rep. 5, 35–54. doi:10.1016/j.bbrep.2015.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Rattu, G., Khansili, N., Maurya, V. K., and Krishna, P. M. (2020). Lactate detection sensors for food, clinical and biological applications: a review. Environ. Chem. Lett. 19, 1135–1152. doi:10.1007/s10311-020-01106-6

CrossRef Full Text | Google Scholar

Rawoof, S. A. A., Kumar, P. S., Vo, D.-V. N., Devaraj, K., Mani, Y., Devaraj, T., et al. (2020). Production of optically pure lactic acid by microbial fermentation: a review. Environ. Chem. Lett. 19, 539–556. doi:10.1007/s10311-020-01083-w

CrossRef Full Text | Google Scholar

Reyné, B., Saby, N., and Sofonea, M. T. (2022). Principles of mathematical epidemiology and compartmental modelling application to covid-19. Anaesth. Crit. Care amp; Pain Med. 41, 101017. doi:10.1016/j.accpm.2021.101017

CrossRef Full Text | Google Scholar

Romero, M. R., Ahumada, F., Garay, F., and Baruzzi, A. M. (2010). Amperometric biosensor for direct blood lactate detection. Anal. Chem. 82, 5568–5572. doi:10.1021/ac1004426

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, S., Sharma, A., and Gupta, M. N. (2006). Preparation of cross-linked enzyme aggregates by using bovine serum albumin as a proteic feeder. Anal. Biochem. 351, 207–213. doi:10.1016/j.ab.2006.01.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Simelevicius, D., Baronas, R., and Kulys, J. (2012). Modelling of amperometric biosensor used for synergistic substrates determination. Sensors 12, 4897–4917. doi:10.3390/s120404897

PubMed Abstract | CrossRef Full Text | Google Scholar

Standards, I. F. (2021). fao.org/fao-who-codexalimentarius/sh-proxy/en/?lnk=1&url=https%253a%252f%252fworkspace. Available at: https://www.fao.org/fao-who-codexalimentarius/sh-proxy/en/?lnk=1&url=https%253A%252F%252Fworkspace.fao.org%252Fsites%252Fcodex%252FStandards%252FCXG%2B36-1989%252FCXG_036e.pdf.

Google Scholar

Suman, S., Singhal, R., Sharma, A. L., Malthotra, B., and Pundir, C. (2005). Development of a lactate biosensor based on conducting copolymer bound lactate oxidase. Sensors Actuators B Chem. 107, 768–772. doi:10.1016/j.snb.2004.12.016

CrossRef Full Text | Google Scholar

Thompson, M. D., and Beard, D. A. (2012). Physiologically based pharmacokinetic tissue compartment model selection in drug development and risk assessment. J. Pharm. Sci. 101, 424–435. doi:10.1002/jps.22768

PubMed Abstract | CrossRef Full Text | Google Scholar

Trifirò, A., Saccani, G., Gherardi, S., Vicini, E., Spotti, E., Previdi, M., et al. (1997). Use of ion chromatography for monitoring microbial spoilage in the fruit juice industry. J. Chromatogr. A 770, 243–252. doi:10.1016/s0021-9673(97)00049-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, Q., and Yeung, E. S. (1994). Indirect fluorescence determination of lactate and pyruvate in single erythrocytes by capillary electrophoresis. J. Chromatogr. A 661, 287–295. doi:10.1016/0021-9673(94)85196-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zouaoui, F., Zine, N., Errachid, A., and Jaffrezic-Renault, N. (2022). Mathematical model and numerical simulation of conductometric biosensor of urea. Electroanalysis 34, 1131–1140. doi:10.1002/elan.202100610

CrossRef Full Text | Google Scholar