Adaptive Estimation of Biological Rhythm in Crassulacean Acid Metabolism with Critical Manifold

The mechanism of endogenous circadian photosynthesis oscillations of plants performing crassulacean acid metabolism (CAM) is investigated in terms of a nonlinear theoretical model. Blasius et al. used throughout continuous time differential equations which adequately reflect the CAM dynamics.Themodel shows regular endogenous limit cycle oscillations that are stable for a wide range of temperatures in a manner that complies well with experimental data. In this paper, we pay attention to the approximation of the fast modes of the CAM dynamics. Using the zero-epsilon approximation of the slow manifold, we derive the critical manifold that is defined by two algebraic nonlinear equations. The critical manifold allows us to give the algebraic estimate of the order of the tonoplast membrane. The dynamic equation of the order of the tonoplast membrane includes the nonlinear function that gives the equilibrium value of the lipid order of tonoplast functions as a hysteresis switch. We identify the nonlinear function with the measurement signals. Using the basis function expansion of the nonlinear and the critical manifold, we propose an adaptive observer to estimate the tonoplast order and the nonlinear function.


Introduction
Biological rhythm is characterized by free-running, endogenous rhythms, ranging from periods of seconds (e.g., heart beat) to years (e.g., population dynamics).Endogenous oscillators of living cells and organisms are based on feedback loops in complex biochemical-reaction networks.Crassulacean acid metabolism (CAM) is a specific mechanism of inorganic carbon acquisition for photosynthesis [1,2].CAM plants take up external CO 2 during the night while the stomata are open.The CO 2 is fixed in the plant in the form of organic acids, normally malate, via the enzyme phosphonenolpyruvate carboxylase (PEPc) and stored overnight in the vacuole [3].The malic acid is accumulated to high concentrations of up to several hundred millimolars and stored in the cell sap vacuole [1,2].During the daytime, malate is remobilized from the vacuole and decarboxylated again [3].The resulting CO 2 is refixed by RubisCO for photosynthetic assimilation [3].Blasius et al. investigated the mechanism of endogenous circadian photosynthesis oscillations of plants performing CAM in terms of a nonlinear theoretical model [3][4][5].They used throughout continuous time differential equations which adequately reflect the CAM dynamics.The model showed regular endogenous oscillations that were stable limit cycles attracting all neighboring trajectories for a wide range of temperatures.
The lack of reliable sensors and the high cost of advanced instrumentation are significant problems in monitoring and control of biological systems.We have recently discussed the nonlinear dynamical model of CAM from the control theoretical viewpoint.The state variables of the nonlinear dynamic equations are an internal CO 2 concentration, , a malate concentration in the cytoplasm, , a malate concentration in the vacuole, , and an order of the tonoplast membrane, , which is called a tonoplast order.The input variables are an external CO 2 concentration, a light intensity, and a temperature.The output is assumed to be a part of the state variables.The order of the tonoplast membrane, , is immeasurable because it is the mean orientational order of lipid chains [6].A software sensor for online estimation of  is needed with available signals.We presented a dynamic estimator of the tonoplast order and a fuzzy identifier of the nonlinear function in the dynamics of the tonoplast order [7] and designed an adaptive observer-based -controller for tonoplast order to synchronize a desired periodic signal under the assumption that the available signals are , , and  [8].The assumption allows us to construct a 2nd-order observer in order to estimate the tonoplast order, , and the nonlinear function in the dynamics of the tonoplast order, (, ) [9].From the viewpoint of real measurements, we presented an adaptive observer to estimate the states and the nonlinear function in the dynamics of the tonoplast order assuming that the available signal is only the internal CO 2 concentration [10].However, the proposed observer with the internal CO 2 concentration does not converge to real states but has other large vibrations in the surrounding of the real oscillation cycles.
In this paper, we design an algebraic estimate of the tonoplast order with the approximation of the critical manifold and a first-order adaptive observer to estimate the tonoplast order and the nonlinear function in the dynamics of the tonoplast order.To begin with, we pay attention to the approximation of the fast modes of the CAM dynamics.Using the zero-epsilon approximation of the slow manifold, that is, the critical manifold, the malate concentration in the cytoplasm, and  is approximated by the nonlinear function of the internal CO 2 concentration, .This means that we can remove the sensor for measuring the malate concentration in the cytoplasm, .Then, we derive an algebraic relationship among the internal CO 2 concentration, , the malate concentration in the vacuole, , and the order of the tonoplast membrane, .This equation is available as an estimate of , when  and  are measurable.We call it the algebraic estimate of .The algebraic estimate does not depend on the temperature.Finally, we design a first-order adaptive observer for online estimation of the nonlinear function in the dynamics of the tonoplast order using the output signal as the algebraic estimate of .The use of the algebraic estimate allows us to improve the convergence speed of the estimator in the previous paper [9,11,12].The observer can attenuate the estimation error of the tonoplast order against the measurement noises.

The Minimal CAM Model
The CAM model that we will use has been studied by Blasius et al. [3,13], and here we only outline the minimal CAM model.The model can be characterized by the block diagram shown in Figure 1 [3,13].The figure shows the major reactant pools of CAM that generate the carbon flow during the circadian cycle.The pool concentrations are the following: (i) internal CO 2 concentration, ; (ii) malate concentration in the cytoplasm, ; (iii) malate concentration in the vacuole, ; (iv)  is a variable that describes the ordering of the lipid molecules in the tonoplast membrane.These are the dynamic variables of the cyclic process.They are connected by the flows,  1 ,  2 , and  3 during the gain and loss terms of the metabolites.The model depends on three external control parameters: temperature, , light intensity, , and external CO 2 concentration,  ext .In this model, regulation of the transport of malate between the cytoplasm and vacuole is the key process for establishing the observed endogenous rhythm in CAM, and it requires a hysteric switching of the passive malate efflux.Blasius et al. [3] introduced  as a new dynamic variable, leading to an additional equation for its time change.The dynamics are characterized by a set of four coupled, nonlinear differential equations of first order in time: where the function (, ) is the thermodynamic equilibrium value of malate concentration in the vacuole, , and is a thirdorder nonlinear function depending on the temperature, , shown in Figure 2 [3].The temperature-dependent -null cline,  = (, ), corresponds to the hysteretic behavior of the phase diagram in the membrane model of Neff et al. [6].The constant, , is the time constant for relaxation into thermal equilibrium.The smallest parameter  ≪ 1 reflects the volume ratio of cytoplasm to vacuole, which is typically of the order of 1/100 in CAM plants [3].The flows   ,  = 1, 2, 3 involve modeling of the metabolic reactions and comprise the whole structure of the carbon circulation in CAM [13].They are described by the following equations: where  1 ,  2 , and  3 are defined as follows [13]: (i)  1 : the difference between malate influx and efflux into and out of the vacuole, modeled with the dynamic hysteresis, (ii)  2 : the difference between malate production from CO 2 fixation by phosphoenolpyruvate carboxylase (PEPc) and its depletion by decarboxylation, (iii)  3 : CO 2 uptake from outside,  ext (), minus CO 2 consumption by photosynthesis, which is directly proportional to the external control parameter light intensity, (), plus CO 2 production by respiration.

Assumptions and Problem Formulation
Throughout this paper, we make the following assumptions.
(i) All states, , , , and , are bounded and positive.
(iii) All parameters are known and the same as in [3].
The order of the tonoplast membrane, , is immeasurable because it is the mean orientational order of lipid chains [6].Moreover, the function form of (, ) has to be numerically derived for each temperature, , since the function (, ) is the thermodynamic equilibrium value of malate concentration in the vacuole [6].However, the computational cost is high [6].Thus, a software sensor for online estimation of  and (, ) is needed with available signals.The design problem is to estimate the tonoplast order, , using the available signals in the case of unknown (, ).
For joint estimation of the states and the nonlinear function, we employ the adaptive observer rather than the Kalman filter.Compared to the Kalman filter, the advantages of the adaptive observer are as follows.
(i) Since the Kalman filter is applied to the extended system, the size of the matrix is larger than that of the adaptive observer.The numerical advantage of the adaptive observer is obvious.
(ii) The stability condition of the adaptive observer is simpler than that of the Kalman filter.

Algebraic Estimate of Tonoplast Order
We can decompose the CAM dynamics with fast and slow variables.The fast subsystem consists of the internal CO 2 concentration, , and the malate concentration in the cytoplasm, .Assuming that  = /, we have the boundary-layer model: The variables  and  in the slow subsystem is injected to the -dynamics as an additional input /.Moreover, the states , , and  do not depend on the temperature  explicitly.To assure the boundedness of the fast subsystem, we assume that () = 1,  ext () = 1, and  and  are constant.
The Jacobian matrix of the boundary-layer model is given by where The trace and the determinant of the Jacobian matrix are calculated as follows: trace Since  1 () is positive for  > 0, the eigenvalues of the Jacobian matrix are negative.Thus, the boundary layer system has an exponentially stable equilibrium point as long as  and  are constant.
In the zero-epsilon approximation ( = 0), the zero-order approximation of the slow manifold, that is, critical manifold [14], is given by The critical manifold is equivalent to the zero-order approximation of the slow manifold [15] and includes the nullcline of -dynamics in the form of the equation: ℎ(, ) = 0, which is independent of the temperature  since  2 and  3 does not depend on .We assume that the equation ℎ(, ) = 0 has the isolated solution  =  1 ().The function  1 (), however, cannot be easily obtained because the expression of  3 is complex.Thus, we get the numerical approximation by MATLAB/Simulink.Figure 4 shows the nullclines of -dynamics and -dynamics.Using  the MATLAB plot editor, we get the following 2nd-order polynomial form as an approximation of the isolated solution  =  1 (): As  depends only on , the value of the malate concentration in the cytoplasm, , can be obtained by the value of the value of the internal CO 2 concentration, , for every temperature.Thus, the function form of the -nullcline is invariant to temperature changes.
When  and  are not constant but oscillating,  and  are oscillating on the equation  =  1 () for large time.In this case, the slow subsystem can be considered as a signal generator for the fast subsystem.Figure 5 shows the phase plots of (, ) for each temperature,  = 0.2246, 0.2254, 0.2242, 0.2238.These lines are on the nullcline (8) except for the neighborhood in the initial value.
Moreover, the right equation of ( 7) is the nullcline of dynamics and is given by Thus, we get the zero-order approximation of the slow manifold as follows: Since the fast subsystem has an exponentially stable equilibrium point, from Tikhonov's theorem [16], we have the following slow manifold: where if  is measurable, an algebraic estimate of  is obtained by using the nonlinear function  1 () as follows: This means that we need not measure  when  is measurable.Moreover,  is not measurable because it is the order parameter of the membrane lipids.If  is measurable,  can be estimated using algebraic equation (11).The algebraic estimate of  is obtained by Defining the estimation error between  and  as , it can be approximated by where  2 is the first-order coefficient of () in the second equation of (11).
Remark 1.If the light intensity, , and the external CO 2 concentration,  ext , are time-varying, -nullcline is given by Since () and  ext () are assumed to be measurable, we can calculate  =  1 () numerically.

Adaptive Observer to Estimate the Nonlinear Function
The nonlinear function (, ) gives the equilibrium value of the lipid order of tonoplast functions as a hysteresis switch that is numerically calculated by the thermodynamics average [6].Since the function form of (, ) has to be numerically derived for each , the computational cost is high.Thus, the soft sensor is needed to estimate the nonlinear function (, ) with the measurement signals.In the limit cycle case such as  = 0.2246, however, we cannot estimate the value of (, ) by using the steady-state value of  and -nullcline because the slow subsystem is not stable and  keeps oscillating.Thus, we need a dynamic estimator with the measurements.
The slow subsystem of the CAM dynamics is given by Since the nonlinear function (, ) is numerically calculated by two integral equations including Boltzman distribution [6], we cannot get the function form of (, ) but obtain a table lookup.We approximate (, ) by the interpolation using the following basis function expansion: where  is the preassigned vector of basis functions,   is the vector of the unknown coefficients for each temperature, and  is the approximation error.This basis function expansion is a kind of the linear parametric form and is widely used to identify nonlinear functions [17,18].In this paper, we adopt the basis functions as the cubic spline interpolation such as (, ) ≈ ∑ 8 =0 ∑ 3 =0   ( −   )  because the spline interpolation is one of the easiest ways to approximate the table lookup.The constants   ,  = 0, . . ., 8 are the interpolation nodes and are selected as (0.028, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4) .
Using the algebraic estimate of , , we obtain the estimate of the nonlinear function (, ) as follows: where θ is the estimate of   obtained by the adaptive observer.Note that  is not measurable because it is the order parameter of the membrane lipids.Using the algebraic estimate , we design the following adaptive observer to estimate (, ): where  is a positive constant.Its adaptive update law is given by the following -modification form: where Γ is a positive definite matrix and  is a positive constant.Denoting the estimation error between ẑ and  as   = ẑ − , we have where θ = θ −   and (, ) = (1/)   (() − ()).Assuming that the function  satisfies the Lipschitz condition we have the following inequality: Defining the Lyapunov function candidate as its derivative satisfies the following inequality:

Simulation Results
In this chapter, we compare the estimation performance of the adaptive observer with the algebraic estimate .Using the parameters as mentioned in Chapter 2, we perform the computer simulation by MATLAB/Simulink when the temperature is  = 0.2246.The initial conditions of the states are given by (0) = 0.4, (0) = 0.62, and (0) = 0.56, (0) = 0.2.The design parameters and the initial condition of the adaptive observer are selected as  = 100, Γ = 50 36 ,  = 0.5, and ẑ(0) = 0.1.Figure 6 shows the algebraic estimate  and the dynamic estimate by the adaptive observer, ẑ, of the tonoplast order.The dynamic estimate achieves almost the same performance.Next, Figure 7 shows the algebraic estimate and the dynamics estimate when the Gaussian white noise is added to the internal CO 2 concentration such as () + V(), where V() ∈ (0, 10 −4 ).The algebraic estimate, , is sensitive to the measurement noise.On the other hand, the dynamic estimate, ẑ, attenuates the vibration as a kind of low-pass filter.Moreover, the adaptive observer allows us to estimate the nonlinear function (, ) as the basis function expansion ĝ(, ) for each temperature.shows the estimate ĝ(, ) in noise-free case, and the right figure shows one in the case where the Gaussian white noise is added to the internal CO 2 concentration such as () + V(), where V() ∈ (0, 10 −4 ).Thus, the adaptive observer is efficient to estimate the value of (, ) with noisy case.The same observer is used to estimate  and (, ) at the other temperatures.The adaptive observer can estimate the dynamics of , and its estimation performance is robust against the temperature changes.

Conclusion
We treated the estimation problem of the CAM plant as a biological system with a limit cycle.The algebraic estimates of the malate concentration in the cytoplasm, , and the tonoplast order  are presented using the critical manifold, which can be used at all temperatures.Moreover, the first-order adaptive observer to estimate the nonlinear function (, ) is designed under the assumption that the signals of the internal CO 2 concentration, , and the malate concentration in the vacuole, , are measurable.The adaptive observer based on the -modification-type update law with the algebraic estimate of  guarantees the boundedness of estimation errors and allows us to reduce the order of the estimator and improve the estimation performance.The adaptive observer can be used as an online soft sensor for the equilibrium value of the lipid order of tonoplast functions as a hysteresis switch.
Using computer simulations, we showed that the estimation performance of the adaptive observer is robust against the measurement noises in .

Figure 3 :
Figure 3: Time responses of sustained endogenous rhythms in continuous light when  = 0.2242.
is negative.Thus, we conclude that   and θ are bounded.The design parameters  and  in the adaptive observer can be used to reduce the estimation errors   and θ .The -modification update law with the algebraic estimate does not need the persistently exciting condition.The proposed adaptive observer is useful in the case where the light intensity, , and the external CO 2 concentration,  ext , are time-varying because the slow subsystem does not depend on  and  ext .