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

Authors: Takami Matsuo; Yusuke Totoki; Haruo Suemitsu
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

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. The model 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.


---

## Body

## 1. 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 CO2 during the night while the stomata are open. The CO2 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 CO2 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–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 CO2 concentration, w, a malate concentration in the cytoplasm, x, a malate concentration in the vacuole, y, and an order of the tonoplast membrane, z, which is called a tonoplast order. The input variables are an external CO2 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, z, is immeasurable because it is the mean orientational order of lipid chains [6]. A software sensor for online estimation of z 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 P-controller for tonoplast order to synchronize a desired periodic signal under the assumption that the available signals are w, x, and y [8]. The assumption allows us to construct a 2nd-order observer in order to estimate the tonoplast order, z, and the nonlinear function in the dynamics of the tonoplast order, g(z,T) [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 CO2 concentration [10]. However, the proposed observer with the internal CO2 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, andx is approximated by the nonlinear function of the internal CO2 concentration, w. This means that we can remove the sensor for measuring the malate concentration in the cytoplasm, x. Then, we derive an algebraic relationship among the internal CO2 concentration, w, the malate concentration in the vacuole, y, and the order of the tonoplast membrane, z. This equation is available as an estimate of z, when w and y are measurable. We call it the algebraic estimate of z. 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 z. 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.


## 2. 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 CO2 concentration, w;(ii)
malate concentration in the cytoplasm,x;(iii)
malate concentration in the vacuole,y;(iv)
z is a variable that describes the ordering of the lipid molecules in the tonoplast membrane.Figure 1
Block diagram of the CAM model showing as dynamic variables (encircled) the three reactant pools,w,  x,  y, and a variable characterizing the order of the lipid molecules in the tonoplast, z (circled), within a CAM cell. They are connected by flows u1, u2, and u3 (solid arrows). Regulatory feedback loops are indicated as dashed lines. The model depends on three external control parameters: temperature, T, light intensity, L, and external CO2 concentration, Cext (dotted arrows) [3].These are the dynamic variables of the cyclic process. They are connected by the flows,u1, u2, and u3 during the gain and loss terms of the metabolites. The model depends on three external control parameters: temperature, T, light intensity, L, and external CO2 concentration, Cext. 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 z 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:
(1)εw˙=-u2+u3,εx˙=-u1+u2,y˙=u1,τz˙=g(z,T)-y,
where the function g(z,T) is the thermodynamic equilibrium value of malate concentration in the vacuole, y, and is a third-order nonlinear function depending on the temperature, T, shown in Figure 2 [3]. The temperature-dependent z-null cline, y=g(z,T), 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].Figure 2
The nonlinear functiony=g(z,T) when T=0.2238, 0.2244, 0.2246, 0.2250, 0.2254.The flowsui, i=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:
(2)u1=cx-yz,u2=wx-x,u3=cJ(Cext(t)-w)exp(αw)-L(t)w+cRLKL(t)+LKw1w+w1,
where u1, u2, and u3 are defined as follows [13]:  (i)
u1: the difference between malate influx and efflux into and out of the vacuole, modeled with the dynamic hysteresis,(ii)
u2: the difference between malate production from CO2 fixation by phosphoenolpyruvate carboxylase (PEPc) and its depletion by decarboxylation,(iii)
u3: CO2 uptake from outside, Cext(t), minus CO2 consumption by photosynthesis, which is directly proportional to the external control parameter light intensity, L(t), plus CO2 production by respiration.Blasius et al. calculated the dynamic behavior using the dimensionless variables with the parameters [3]: Cext=1, L(t)=1, T=0.2238, 0.2242, 0.2246, 0.2250, 0.2254, c=5.5, cJ=1, cR=1, ε=0.001, τ=0.35, α=1.5, w1=0.1, LK=0.5, and R=0.1. Using these parameters, we perform the computer simulation by MATLAB/Simulink. The initial conditions of the states are given as w(0)=0.4, x(0)=0.62, y(0)=0.56, and z(0)=0.2. The nonlinear function g(z,T) is approximated via the third-order spline interpolation from the figure shown in [3]. Figure 3 shows the dynamic responses of the state variables when T=0.2242. Blasius et al. have shown the influence of temperature as an external control parameter [3]. With increasing temperature, the steady state malate level decreases, until, at the critical temperature, the system undergoes supercritical Hopf bifurcation, and the stable fixed point changes into an unstable one surrounded by a small-amplitude limit cycle [3].Figure 3
Time responses of sustained endogenous rhythms in continuous light whenT=0.2242.


## 3. Assumptions and Problem Formulation

Throughout this paper, we make the following assumptions.(i)
All states,w, x, y, and z, are bounded and positive.(ii)
The available signals arew, y, T, L, and Cext.(iii)
All parameters are known and the same as in [3].The order of the tonoplast membrane,z, is immeasurable because it is the mean orientational order of lipid chains [6]. Moreover, the function form of g(z,T) has to be numerically derived for each temperature, T, since the function g(z,T) 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 z and g(z,T) is needed with available signals. The design problem is to estimate the tonoplast order, z, using the available signals in the case of unknown g(z,T).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.


## 4. Algebraic Estimate of Tonoplast Order

We can decompose the CAM dynamics with fast and slow variables. The fast subsystem consists of the internal CO2 concentration, w, and the malate concentration in the cytoplasm, x. Assuming that l=t/ε, we have the boundary-layer model:
(3)dwdl=-u2+u3,dxdl=-u1+u2.The variablesy and z in the slow subsystem is injected to the x-dynamics as an additional input y/z. Moreover, the states w, x, and y do not depend on the temperature T explicitly. To assure the boundedness of the fast subsystem, we assume that L(t)=1, Cext(t)=1, and y and z are constant.The Jacobian matrix of the boundary-layer model is given by(4)Jwx=[J11J12J21J22]=[∂(-u2+u3)∂w∂(-u2+u3)∂x∂(-u1+u2)∂w∂(-u1+u2)∂x],J11=-1x-a1(w),J12=wx2+1,J21=1x,J22=-wx2-(c+1),
where
(5)a1(w)=cJ1+(1-w)αexp(αw)+1+cRLK1+LKw1(w+w1)2.The trace and the determinant of the Jacobian matrix are calculated as follows:(6)traceJwx=-1x-a1(w)-wx2-(c+1),detJwx=a1(w)(1+wx2)+c(a1(w)+1x).Sincea1(w) is positive for w>0, the eigenvalues of the Jacobian matrix are negative. Thus, the boundary layer system has an exponentially stable equilibrium point as long as y and z are constant.In the zero-epsilon approximation (ε=0), the zero-order approximation of the slow manifold, that is, critical manifold [14], is given by
(7)0=-u2+u3,0=-u1+u2.The critical manifold is equivalent to the zero-order approximation of the slow manifold [15] and includes the nullcline of w-dynamics in the form of the equation: h(w,x)=0, which is independent of the temperature T since u2 and u3 does not depend on T. We assume that the equation h(w-,x-)=0 has the isolated solution x-=ϕ1(w-). The function ϕ1(w-), however, cannot be easily obtained because the expression of u3 is complex. Thus, we get the numerical approximation by MATLAB/Simulink. Figure 4 shows the nullclines of w-dynamics and x-dynamics. Using the MATLAB plot editor, we get the following 2nd-order polynomial form as an approximation of the isolated solution x-=ϕ1(w-):
(8)x-≈-0.43146w-2+2.2032w--0.17466=:ϕ1(w-).Figure 4
w-nullcline and x-nullcline.Asx- depends only on w-, the value of the malate concentration in the cytoplasm, x, can be obtained by the value of the value of the internal CO2 concentration, w, for every temperature. Thus, the function form of the w-nullcline is invariant to temperature changes.Wheny and z are not constant but oscillating, x and w are oscillating on the equation x=ϕ1(w) 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 (w,z) for each temperature, T=0.2246,0.2254,0.2242,0.2238. These lines are on the w-nullcline (8) except for the neighborhood in the initial value.The phase plane(w,z) ((a) T=0.2246, (b) T=0.2254, (c) T=0.2242, and (d) T=0.2238).
(a)
(b)
(c)
(d)Moreover, the right equation of (7) is the nullcline of x-dynamics and is given by
(9)yz=(1+c)x--w-x-=:ϕ2(w-,x-).Thus, we get the zero-order approximation of the slow manifold as follows:(10)x-=ϕ1(w-),yz=ϕ2(w-,x-).Since the fast subsystem has an exponentially stable equilibrium point, from Tikhonov's theorem [16], we have the following slow manifold:
(11)x=ϕ1(w)+O(ε),yz=ϕ2(w,x)+O(ε),
where if w is measurable, an algebraic estimate of x is obtained by using the nonlinear function ϕ1(w) as follows:
(12)x≈ϕ1(w).This means that we need not measurex when w is measurable.Moreover,z is not measurable because it is the order parameter of the membrane lipids. If y is measurable, z can be estimated using algebraic equation (11). The algebraic estimate of z is obtained by
(13)z-=yϕ2(w,ϕ1(w)).Defining the estimation error betweenz and z- as ϵ, it can be approximated by
(14)ϵ=z--z≈yϕ2-yϕ2+L2ε≈εyL2ϕ22,
where L2 is the first-order coefficient of O(ε) in the second equation of (11).Remark 1.
If the light intensity,L, and the external CO2 concentration, Cext, are time-varying, w-nullcline is given by
(15)-wx+x+cJ(Cext(t)-w)exp(αw)-L(t)w+cRLKL(t)+LKw1w+w1=0.
Since L(t) and Cext(t) are assumed to be measurable, we can calculate x-=ϕ1(w-) numerically.


## 5. Adaptive Observer to Estimate the Nonlinear Function

The nonlinear functiong(z,T) 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 g(z,T) has to be numerically derived for each T, the computational cost is high. Thus, the soft sensor is needed to estimate the nonlinear function g(z,T) with the measurement signals. In the limit cycle case such as T=0.2246, however, we cannot estimate the value of g(z,T) by using the steady-state value of y and z-nullcline because the slow subsystem is not stable and z keeps oscillating. Thus, we need a dynamic estimator with the measurements.The slow subsystem of the CAM dynamics is given by(16)y˙=cx-yz,τz˙=g(z,T)-y.Since the nonlinear functiong(z,T) is numerically calculated by two integral equations including Boltzman distribution [6], we cannot get the function form of g(z,T) but obtain a table lookup. We approximate g(z,T) by the interpolation using the following basis function expansion:
(17)g(z,T)=θT′ζ(z)+δ,
where ζ is the preassigned vector of basis functions, θT 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 g(z,T)≈∑i=08‍∑j=03‍θij(z-ai)j because the spline interpolation is one of the easiest ways to approximate the table lookup. The constants ai, i=0,…,8 are the interpolation nodes and are selected as
(18)(0.028,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4).Using the algebraic estimate ofz, z-, we obtain the estimate of the nonlinear function g(z,T) as follows:
(19)g^(z-,T)=θ^T′(t)ζ(z-),
where θ^T is the estimate of θT obtained by the adaptive observer.Note thatz is not measurable because it is the order parameter of the membrane lipids. Using the algebraic estimate z-, we design the following adaptive observer to estimate g(z,T):
(20)z^˙=1τ{θ^T′ζ(z-)-y}+k(z--z^),
where k is a positive constant. Its adaptive update law is given by the following σ-modification form:(21)θ^˙T=-Γζ(z-)(z^-z-)-σθ^,
where Γ is a positive definite matrix and α is a positive constant. Denoting the estimation error between z^ and z as ez=z^-z, we have
(22)e˙z=-kez+kϵ+1τθ~T′ζ(z-)+φ(z-,z,T)-δτ,
where θ~T=θ^T-θT and φ(z-,z)=(1/τ)θT′(ζ(z-)-ζ(z)). Assuming that the function ζ satisfies the Lipschitz condition
(23)∥ζ(z-)-ζ(z)∥≤Kζ|ϵ|,
we have the following inequality:
(24)|φ(z-,z)|≤1τKζ∥θT∥|ϵ|.Defining the Lyapunov function candidate as(25)V=12(ez2+1τθ~T′Γ-1θ~T),
its derivative satisfies the following inequality:
(26)V˙=-kez2-στθ~T′θ~T-1τezθT′(ζ(z-)-ζ(z))+1τθ~T′ζ(z-)ϵ-δτez≤-kez2-στθ~T′θ~T+1τ∥θT∥|ez|Kζ|ϵ|+στ|ez|+1τ∥θ~T∥|ζ(z-)||ϵ|=-|ez|(k|ez|-δ+Kζ∥θT∥|ϵ|τ)-στ∥θ~T∥(∥θ~T∥-1σ|ζ(z-)||ϵ|).Note thatϵ does not depend on the adaptive observer but depends on the algebraic estimate. If the following equations are satisfied:
(27)|ez|>δ+Kζ∥θT∥|ϵ|kτ,∥θ~T∥>1σ|ζ(z-)||ϵ|,
then V˙ is negative. Thus, we conclude that ez and θ~ are bounded. The design parameters k and σ in the adaptive observer can be used to reduce the estimation errors ez and θ~T. The σ-modification update law with the algebraic estimate does not need the persistently exciting condition.Remark 2.
The proposed adaptive observer is useful in the case where the light intensity,L, and the external CO2 concentration, Cext, are time-varying because the slow subsystem does not depend on L and Cext.


## 6. Simulation Results

In this chapter, we compare the estimation performance of the adaptive observer with the algebraic estimatez-. Using the parameters as mentioned in Chapter 2, we perform the computer simulation by MATLAB/Simulink when the temperature is T=0.2246. The initial conditions of the states are given by w(0)=0.4, x(0)=0.62, and y(0)=0.56, z(0)=0.2. The design parameters and the initial condition of the adaptive observer are selected as k=100, Γ=50I36, σ=0.5, and z^(0)=0.1. Figure 6 shows the algebraic estimate z- and the dynamic estimate by the adaptive observer, z^, 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 CO2 concentration such as w(t)+v(t), where v(t)∈N(0,10-4). The algebraic estimate, z-, is sensitive to the measurement noise. On the other hand, the dynamic estimate, z^, attenuates the vibration as a kind of low-pass filter. Moreover, the adaptive observer allows us to estimate the nonlinear function g(z,T) as the basis function expansion g^(z-,T) for each temperature. Figure 8 shows g(z,T) and g^(z-,T) when T=0.2246. The left figure shows the estimate g^(z-,T) in noise-free case, and the right figure shows one in the case where the Gaussian white noise is added to the internal CO2 concentration such as w(t)+v(t), where v(t)∈N(0,10-4). Thus, the adaptive observer is efficient to estimate the value of g(z,T) with noisy case. The same observer is used to estimate z and g(z,T) at the other temperatures. The adaptive observer can estimate the dynamics of z, and its estimation performance is robust against the temperature changes.The responses ofz and its estimates z-, z^ when T=0.2246 ((a) z (solid line) and z- (dotted line), (b) z (solid line) and z^ (dotted line)).
(a)
(b)The responses ofz and its estimates z-, z^ when T=0.2246 ((a) z (solid line) and z- (dotted line), (b) z (solid line) and z^ (dotted line)).
(a)
(b)The estimated nonlinear functiong^(z-,T) when T=0.2246.
(a)
(b)


## 7. 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,x, and the tonoplast order z are presented using the critical manifold, which can be used at all temperatures. Moreover, the first-order adaptive observer to estimate the nonlinear function g(z,T) is designed under the assumption that the signals of the internal CO2 concentration, w, and the malate concentration in the vacuole, y, are measurable. The adaptive observer based on the σ-modification-type update law with the algebraic estimate of z 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 w.



---

*Source: 856404-2013-06-18.xml*