Abstract
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” ().
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) () 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 ().
FIGURE 1
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 (). In addition, lactate can act as an indicator of the activity of bacteria during the fermentation process (). It can indicate the freshness and quality of some products—wine (), juices (), 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 (; ), spectrophotometry (), fluorometry (), liquid chromatography (), and nuclear magnetic resonance (; ). 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 ().
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) (). 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 (). 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 ().
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 ), 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 () and electrochemical—amperometric, potentiometric, or, less often, conductometric or impedimetric (). Moreover, various methods of improving the main characteristics of sensors have been applied to known biosensor systems—nanoparticles (; ), other nanomaterials (), complex, multi-stage methods of enzyme immobilization (), and multi-enzyme membranes.
Although amperometric () and potentiometric () biosensors are usually monoenzymatic, that is, based on LHD or Lx, conductometric () and impedimetric () biosensors usually use two-enzyme bioselective elements, for example, based on a mixture of lactate oxidase and peroxidase () or a mixture of lactate dehydrogenase and pyruvate oxidase ().
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 k1. The resulting complex can dissociate into E and S (with a rate constant of k−1) or transform into a product P with a rate constant of k2 ().
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):Here, v is the rate of reaction (mol/s), Vmax is the maximum reaction rate (mol/s), nS is substrate concentration (mol/dm3), and KM is the Michaelis–Menten constant (mol/dm3).
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 (; ).
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 (). 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 (; ), mathematical epidemiology (; ; ), enzyme kinetics (; ; ), demography (), and ecology (). Stability analysis plays a pivotal role in understanding the behavior of these systems under different conditions (). In recent years, researchers have made significant strides in advancing stability research methodologies, particularly in the context of compartmental systems with delay (; ).
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 (). 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 ().
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
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
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 ), where complex C has a lifetime τ before being decayed. We called the reaction schemean irreversible one-complex Brown’s (IR1CB) mechanism. In ), we offered the following model based on continuously distributed delays:where for confidence level c ∈ (0, 1), we set , f(s) is the density function of the delay distribution, which was designed in the form of a gamma distribution: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 . 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 ) 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 () (Section 1.4)
by adding the time durations τ1 and τ2 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:which we call the irreversible one-complex with two delays Michaelis–Menthen (IR1C2DMM) mechanism. Based on the general approach described by , it yields the delayed model
For the solutions of Eq. 4, elements of which are the vector functions nS, nE, nC, nP ∈ C1 ([−τmax], 0], R4), we consider the following initial conditions:
The value of nP(t) can be found by direct integration, namely,
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 nP(t) as
In , 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 . 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 and .
In the following, we will denote product concentrations obtained as the solutions of the models as nP,pred(t). In turn, the corresponding values of the current due to Eq. 7 will be Ipred(t). On the other hand, let Iexp(t) be the values of responses received as a result of the experiments.
The proposed parameter identification uses the currents Iexp,j (ti) and Ipred,j (ti), , experimentally and numerically obtained at the time instances ti, 1 for given initial substrate concentrations nS (0) = nS,j.
Our goal is to estimate the parameters , when applying Model (2) or (4), respectively.
We will estimate the parameters Π of the models solving the problems of optimization:
Here,is the target function, andare 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 ). 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 τ1 and τ2 when estimating Eq. 4. Note that the mean value of time delay for Model (2) can be calculated as
So, we set the initial values of delays from Eq. 4 such that .
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 formallowing us to approximate the characteristic quasi-polynomial with the help of rational functions ().
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 .
To begin, we gain a thorough understanding of Model (4), including its equations and parameters. We select parameters τ1 and τ2 to vary and the variables nS, nE, and nC 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 nE = d in the phase space, where d = (mintnE(t) + maxtnE(t))/2. This section is determined by specific criteria; namely, we choose such points such that crossing a plane nE = d will happen at sequential time instances t⋆ + T, t⋆ + 2T, t⋆ + 3T, … , where T is a period value; that is, . We plot the sampled points in an (nS, nE)-space to visualize the Poincaré sections.
Repeating this process for various parameter values (τ1, τ2), 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 nS, nE and (nS, nP) 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 nS(0) equal to 0.1 mM, 0.5 mM, 1.0 mM, and 2.5 mM sequentially.
The goal is to estimate the parameters . 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
| kd | a | m | τmin | K | ||
|---|---|---|---|---|---|---|
| Πlower | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 |
| Πupper | 1 | 1,000 | 1,000 | 1,000 | 1.e6 | 1.e6 |
| Πinit | 0.1 | 1 | 20 | 5 | 1.e−10 | 0.5 |
| Πopt | 0.09647381 | 0.9180534 | 16.2354 | 4.607475 | 1e−10 | 0.06818154 |
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
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
3.2 Parameter identification for a model with two discrete delays
Model (4) requires the estimating the parameters . 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
| k1 | k−1 | k2 | τ1 | K | τ2 | ||
|---|---|---|---|---|---|---|---|
| Πlower | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 | 1.e−10 | 0 |
| Πupper | 1 | 1,000 | 1,000 | 1,000 | 1.e6 | 1.e6 | 1,000 |
| Πinit | 0.1 | 0.1 | 0.1 | 15 | 0.5 | 0.0 | 7 |
| Πopt | 0.08717502 | 0.1048326 | 0.1220813 | 15.69849 | 1.171464e-05 | 0.2323843 | 5.357829 |
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
3.3 Qualitative analysis
3.3.1 Existence and positiveness of the solutions
Given , because the right-hand sides of Eq. 4 imply the Lipschitz condition, there exists a unique trajectory of Eq. 4 starting from Eq. 5 ().
Henceforth, we will focus our attention on the positiveness of the solution of the system shown in Eq. 4.
The positiveness of nP(t) follows directly from the positiveness of nC(t). So, we prove the positiveness of nS(t), nE(t), and nC(t) by contradiction.
3.3.1.1 Case without delays
First, we demonstrate the positiveness for Model (4) without delays, that is, if τ1 = τ2 = 0
Let us assume, for the sake of contradiction, that there is the smallest value among tc, te, and ts delivering non-positive solutions. Consider them sequentially.
Let tc > 0 be the smallest instance of time such that nC(tc) = 0. From the first and second lines of Eq. 12, we getIt implies thatFrom the third line of Eq. 12, we haveHence,By continuing Eq. 13, we have nC(tc) > 0, which contradicts the initial assumption.
Let te be the smallest instant that nE(te) = 0. From the third part of Eq. 12, it follows thatand nC(t) > 0 for t ∈ [0, te).
In turn, from the second part of Eq. 12, we haveandfor t ∈ [0, te],
Let ts be the smallest instant that nS (ts) = 0. From the second line of Eq. 12, we get that nC(t) > 0, t ∈ [0, ts). Furthermore, from the first equation, we haveandfor t ∈ [0, ts], 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 (nS(0), nE(0), nC(0), nP(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 τ1 and τ2. 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 ()where Θ(⋅) is the Heaviside function. This function has the property that . Note that contrary to the “undelayed” exponential function, can also accept negative values.
To apply the approach of contradictions shown above, let tc, te, and ts be instances delivering non-positive solutions.
If tc > 0 be the smallest instance that nC (tc) = 0, then from the first and second lines of Eq. 4
Hence,where .
From the third line of Eq. 4, we have
Consider t = tc. Hence,and we see that
Following the proof in the previous case, we conclude that the solutions of Eq. 4 are positive if rate parameters, initial conditions, and delays τ1, τ2 are such that for any t > 0
3.3.2 Stability research
3.3.2.1 Equilibrium of the system
Let be the equilibrium of Model (4). It should satisfyIt follows that
When analyzing Eq. 4, we see that there is a conserved quantity of enzyme because . Hence, , where is the total amount of available enzyme. Thus, from the last equality of Section 3.3.2.1, we conclude that and .
So, we have the unique equilibrium of Model (4), which is the substrate and complex free equilibrium (SCFE).
The value of is bounded and can be determined by Eq. 6. Note that it is related to the initial values of both nS and nE. nP 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
The characteristic quasi-polynomial for is given bywhereHence, stability analysis can be reduced to obtaining conditions of non-positive values of real parts of the quasi-polynomial of χ1(λ).
For SCFE, we have the following quasi-polynomial:
When applying the Padé approximation shown in Eq. 11 to the characteristic quasi-polynomial, we get
Setting the values for , k1, k2, k−1 from Πopt (Table 2), we obtain
In turn, applying the Padé approximation yields
Consider two special cases of Eq. 4:
Case 1 without delays, that is, τ1 = τ2 = 0. Then,
and we have the roots of χ1(λ):
which means marginal stability in this case.
Case 2 with optimal values of the delays from Πopt (Table 2), that is, τ1 = 15.69849 and τ2 = 5.357829. It holds
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 τ1 on Model 4’s dynamics, provided that τ2 = 5.357829 is fixed. To explore this influence, we varied the parameter τ1 over the interval [0, 35].
Figure 7 illustrates the construction of Poincaré sections in 3D plots for various values of τ1. 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 τ1 and τ2 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
Furthermore, Figure 8 presents the corresponding Poincaré sections plotted in the (nS, nE)-plane, offering a clearer visualization of the system’s trajectory crossings. A thorough inspection of the sections says that for the initial trajectory (nS(t), nE(t), nC(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
In Figure 9, we depict the bifurcation diagram for the variable nS resulting from the parameter variation of τ1. It reveals that the equilibrium state behaves as a stable node within the range of τ1 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 (nS, nE)-phase plane for different initial states.
FIGURE 9
FIGURE 10
For τ1 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 nS (0) lead the solution nP to converge toward either a stable node or a limited cycle.
FIGURE 11
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 τ1 + τ2 ≈ E(τ), where τ1, τ2 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 (k1, k2, k−1), time delays τ1, τ2, and the initial concentration of .
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, , nC ≡ 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 nC(t), that is, . In addition, is determined by the initial value nS (0).
In the work, we also considered a model with two discrete delays where τ1 and τ2 are time durations required for the entire fulfillment of forwarded reactions. To be precise, τ1 is the time needed for the substrate to bind with the enzyme (time needed for complex formation), and τ2 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 and and they correspond to a stable solution. If , more time is required for forming the complex at . In other words, less complex will be produced compared with the optimal value. In turn, if , this means that more time will be needed for the 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 τ2 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 and are not constant. They likely depend on the concentrations of S, E (for τ1), and C (for τ2). 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 (). 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.
Statements
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
1.^Without loss of generality, we assume that all time series have the same time span.
References
1
AngstmannC. N.BurneyS.-J. M.HenryB. I.JacobsB. A.XuZ. (2023). A systematic approach to delay functions. Mathematics11, 4526. 10.3390/math11214526
2
AzzouziS.RotariuL.BenitoA. M.MaserW. K.AliM. B.BalaC. (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. 10.1016/j.bios.2015.03.012
3
BakerG. A.Graves-MorrisP. (1996). Padé approximants. Second Edition. Cambridge: Cambridge University Press. 10.1017/cbo9780511530074
4
BergJ.TymoczkoJ.StryerL. (2002). Biochemistry. Fifth Edition. Germany: W. H. Freeman.
5
BiagiS.GhimentiS.OnorM.BramantiE. (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. 10.1002/bmc.2713
6
BlanchiniF.BolzernP.ColaneriP.De NicolaoG.GiordanoG. (2023). Optimal control of compartmental models: the exact solution. Automatica147, 110680. 10.1016/j.automatica.2022.110680
7
BrownA. J. (1902). Xxxvi.—enzyme action. J. Chem. Soc. Trans.81, 373–388. 10.1039/ct9028100373
8
ChanD.BarsanM. M.KorpanY.BrettC. M. (2017). L-lactate selective impedimetric bienzymatic biosensor based on lactate dehydrogenase and pyruvate oxidase. Electrochimica Acta231, 209–215. 10.1016/j.electacta.2017.02.050
9
ChatterjeeY.BhowalB.GuptaK. J.PareekA.Singla-PareekS. L. (2023). Lactate dehydrogenase superfamily in rice and arabidopsis: understanding the molecular evolution and structural diversity. Int. J. Mol. Sci.24, 5900. 10.3390/ijms24065900
10
ChaubeyA.PandeK.SinghV.MalhotraB. (2000). Co-immobilization of lactate oxidase and lactate dehydrogenase on conducting polyaniline films. Anal. Chim. Acta407, 97–103. 10.1016/s0003-2670(99)00797-7
11
CraciunG.MinchevaM.PanteaC.YuP. Y. (2020). Delay stability of reaction systems. Math. Biosci.326, 108387. 10.1016/j.mbs.2020.108387
12
CuiX.LiC. M.ZangJ.YuS. (2007). Highly sensitive lactate biosensor by engineering chitosan/PVI-os/CNT/LOD network nanocomposite. Biosens. Bioelectron.22, 3288–3292. 10.1016/j.bios.2007.03.004
13
GamellaM.CampuzanoS.ConzueloF.CurielJ.MuñozR.ReviejoA.et al (2010). Integrated multienzyme electrochemical biosensors for monitoring malolactic fermentation in wines. Talanta81, 925–933. 10.1016/j.talanta.2010.01.038
14
GibsonT. D. (1999). Biosensors: the stabilité problem. Analusis27, 630–638. 10.1051/analusis:1999270630
15
HaleJ. K.LunelS. M. V. (2013). Introduction to functional differential equations, 99. Germany: Springer Science and Business Media.
16
Hashem ZadehS. M.HeidarshenasM.GhalambazM.NoghrehabadiA.Saffari PourM. (2020). Numerical modeling and investigation of amperometric biosensors with perforated membranes. Sensors20, 2910. 10.3390/s20102910
17
HuS.YoshiharaH. A.BokR.ZhouJ.ZhuM.KurhanewiczJ.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. Imaging30, 1367–1372. 10.1016/j.mri.2012.05.012
18
KeenerJ.SneydJ. (2009). Mathematical physiology (New York: Springer). 10.1007/978-0-387-75847-3
19
KrishnaS.PetersonV.ListmannL.HinnersJ. (2024). Interactive effects of viral lysis and warming in a coastal ocean identified from an idealized ecosystem model. Ecol. Model.487, 110550. 10.1016/j.ecolmodel.2023.110550
20
LockridgeO.MasseyV.SullivanP. A. (1972). Mechanism of action of the flavoenzyme lactate oxidase. J. Biol. Chem.247, 8097–8106. 10.1016/s0021-9258(20)81814-6
21
MartsenyukV.AugustynekK.UrbasA. (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. 10.1016/j.ijnonlinmec.2020.103630
22
MartsenyukV.BernasM.Klos-WitkowskaA. (2021b). Two-strain covid-19 model using delayed dynamic system and big data. IEEE Access9, 113866–113878. 10.1109/access.2021.3104519
23
MartsenyukV.Klos-WitkowskaA.DzyadevychS.SverstiukA. (2022). Nonlinear analytics for electrochemical biosensor design using enzyme aggregates and delayed mass action. Sensors22, 980. 10.3390/s22030980
24
MartsenyukV. P.AndrushchakI. E.MelenchukI. B. (2012). Method of construction and determination of approximate solutions of the model of pharmacokinetics of nanoparticles. J. Automation Inf. Sci.44, 32–43. 10.1615/jautomatinfscien.v44.i8.40
25
MartsenyukV. P.AndrushchakI. Y.GandzyukN. M. (2013). Constructing exponential estimates in compartmental systems with distributed delays: an approach based on the hale–lunel inequality. Cybern. Syst. Analysis49, 347–352. 10.1007/s10559-013-9517-0
26
MartsenyukV. P.GandzyukN. M. (2013). Stability estimation method for compartmental models with delay. Cybern. Syst. Analysis49, 81–85. 10.1007/s10559-013-9488-1
27
MengardaP.DiasF. A.PeixotoJ. V.OsieckiR.BergaminiM. F.Marcolino-JuniorL. 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. 10.1016/j.snb.2019.126663
28
Navarro ValenciaV. A.DíazY.PascaleJ. M.BoniM. F.Sanchez-GalanJ. 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. Heliyon9, e15424. 10.1016/j.heliyon.2023.e15424
29
NesakumarN.SethuramanS.KrishnanU. M.RayappanJ. B. B. (2013). Fabrication of lactate biosensor based on lactate dehydrogenase immobilized on cerium oxide nanoparticles. J. Colloid Interface Sci.410, 158–164. 10.1016/j.jcis.2013.08.009
30
Nguyen-BoisseT. T.SaulnierJ.Jaffrezic-RenaultN.LagardeF. (2013). Highly sensitive conductometric biosensors for total lactate, d- and l-lactate determination in dairy products. Sensors Actuators B Chem.179, 232–239. 10.1016/j.snb.2012.10.021
31
ParkJ. M.JosanS.MayerD.HurdR. E.ChungY.BendahanD.et al (2015). Hyperpolarized 13c NMR observation of lactate kinetics in skeletal muscle. J. Exp. Biol.218, 3308–3318. 10.1242/jeb.123141
32
ParraA.CaseroE.VázquezL.ParienteF.LorenzoE. (2006). Design and characterization of a lactate biosensor based on immobilized lactate oxidase onto gold surfaces. Anal. Chim. Acta555, 308–315. 10.1016/j.aca.2005.09.025
33
PundirC. S.NarwalV. (2017). A bird’s eye view of lactate biosensors. J. Intensive Crit. Care03. 10.21767/2471-8505.100089
34
PundirC. S.NarwalV.BatraB. (2016). Determination of lactic acid with special emphasis on biosensing methods: a review. Biosens. Bioelectron.86, 777–790. 10.1016/j.bios.2016.07.076
35
RabinowitzJ. D.EnerbäckS. (2020). Lactate: the ugly duckling of energy metabolism. Nat. Metab.2, 566–571. 10.1038/s42255-020-0243-4
36
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).
37
RassaeiL.OlthuisW.TsujimuraS.SudhölterE. J. R.van den BergA. (2013). Lactate biosensors: current status and outlook. Anal. Bioanal. Chem.406, 123–137. 10.1007/s00216-013-7307-1
38
RatheeK.DhullV.DhullR.SinghS. (2016). Biosensors based on electrochemical lactate detection: a comprehensive review. Biochem. Biophysics Rep.5, 35–54. 10.1016/j.bbrep.2015.11.010
39
RattuG.KhansiliN.MauryaV. K.KrishnaP. M. (2020). Lactate detection sensors for food, clinical and biological applications: a review. Environ. Chem. Lett.19, 1135–1152. 10.1007/s10311-020-01106-6
40
RawoofS. A. A.KumarP. S.VoD.-V. N.DevarajK.ManiY.DevarajT.et al (2020). Production of optically pure lactic acid by microbial fermentation: a review. Environ. Chem. Lett.19, 539–556. 10.1007/s10311-020-01083-w
41
ReynéB.SabyN.SofoneaM. T. (2022). Principles of mathematical epidemiology and compartmental modelling application to covid-19. Anaesth. Crit. Care amp; Pain Med.41, 101017. 10.1016/j.accpm.2021.101017
42
RomeroM. R.AhumadaF.GarayF.BaruzziA. M. (2010). Amperometric biosensor for direct blood lactate detection. Anal. Chem.82, 5568–5572. 10.1021/ac1004426
43
ShahS.SharmaA.GuptaM. N. (2006). Preparation of cross-linked enzyme aggregates by using bovine serum albumin as a proteic feeder. Anal. Biochem.351, 207–213. 10.1016/j.ab.2006.01.028
44
SimeleviciusD.BaronasR.KulysJ. (2012). Modelling of amperometric biosensor used for synergistic substrates determination. Sensors12, 4897–4917. 10.3390/s120404897
45
StandardsI. 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.
46
SumanS.SinghalR.SharmaA. L.MalthotraB.PundirC. (2005). Development of a lactate biosensor based on conducting copolymer bound lactate oxidase. Sensors Actuators B Chem.107, 768–772. 10.1016/j.snb.2004.12.016
47
ThompsonM. D.BeardD. A. (2012). Physiologically based pharmacokinetic tissue compartment model selection in drug development and risk assessment. J. Pharm. Sci.101, 424–435. 10.1002/jps.22768
48
TrifiròA.SaccaniG.GherardiS.ViciniE.SpottiE.PrevidiM.et al (1997). Use of ion chromatography for monitoring microbial spoilage in the fruit juice industry. J. Chromatogr. A770, 243–252. 10.1016/s0021-9673(97)00049-6
49
XueQ.YeungE. S. (1994). Indirect fluorescence determination of lactate and pyruvate in single erythrocytes by capillary electrophoresis. J. Chromatogr. A661, 287–295. 10.1016/0021-9673(94)85196-4
50
ZouaouiF.ZineN.ErrachidA.Jaffrezic-RenaultN. (2022). Mathematical model and numerical simulation of conductometric biosensor of urea. Electroanalysis34, 1131–1140. 10.1002/elan.202100610
Summary
Keywords
lactate biosensor, operational stability, Michaelis–Menten kinetics, delay, quasi-polynomial, marginal stability, Poincaré section, bifurcation diagram
Citation
Martsenyuk V, Soldatkin O, Klos-Witkowska A, Sverstiuk A and Berketa K (2024) Operational stability study of lactate biosensors: modeling, parameter identification, and stability analysis. Front. Bioeng. Biotechnol. 12:1385459. doi: 10.3389/fbioe.2024.1385459
Received
12 February 2024
Accepted
03 June 2024
Published
16 July 2024
Volume
12 - 2024
Edited by
Kang Cui, University of Jinan, China
Reviewed by
Hui Wang, University of Science and Technology Beijing, China
Lingyin Meng, Linköping University, Sweden
Updates
Copyright
© 2024 Martsenyuk, Soldatkin, Klos-Witkowska, Sverstiuk and Berketa .
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Vasyl Martsenyuk , vmartsenyuk@ubb.edu.pl
Disclaimer
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.