3.1 Cyclic Feedback Systems

A set of n coupled ordinary-differential-equations:

(1)

where:

is called a n–th order cyclic feedback system (see for instance [¹¹] and the references therein) if there exist real constants δ_i ∈ {±1}, i = 1, 2, . . . , n, such that:

Remark 1: As is obvious from the description (1), the state variable x_i−1 drives the state variable x_i. This regulatory characteristic fixes the cyclic feedback defining property of the system.

Cyclic biochemical reactions

Figure 1 shows the particular case of cyclic feed-back systems described by (see for instance [²]):

(2)

Fig. 1 A particular class of cyclic feedback systems. In this case, the rate of change of each state in the chain results from the addition of a function which depends on the state itself plus a function which depends on the contiguous downstream state 

where a_i (·) and b_i (·) are continuous functions (if these functions are real constants the corresponding system is obviously linear).

Remark 2: The cyclic feedback modular structure shown in (2) has been frequently considered to describe cyclic biochemical reactions, where the end product drives the first reaction, and with the functions a_i (·) and b_i (·) usually related to enzymatic processes. The cyclic feedback modular structure (2) is particularly suited for both the description and the design of biosynthetic oscillatory systems. Notice that for biochemical reactions each state takes positive values, i.e., x_i ∈ ℜ⁺.

3.2 Synthesized Cyclic Feedback Transcriptional Networks

Biosynthetic oscillators based on transcriptional regulation are frequently built following the cyclic structure given by (2), with inhibition of gene expression as the basic regulatory mechanism (see for instance [¹³]). We must point out that common natural cell biochemical oscillators (e.g. circadian oscillators in animals and fungus), involve both transcriptional regulation and translational regulation (see for instance [¹⁴]).

Since the cell transcriptional machinery is easy to manipulate via standard recombinant DNA technologies, feedback-based transcriptional oscillators are preferred in cell biosynthetical constructs. We describe here the particular class of cyclic feedback systems that concerns our proposal. We focus our attention in cyclic feedback biosynthetic networks resulting from the chaining of transcriptional modules. Therefore, we first introduce the basic representation of a transcriptional module in input-output terms.

Consider x to be a gene, and let us assume that it codes the transcription factor¹ [x], where [∗] stands for the concentration of the protein coded by ∗. We can now recall the following (see for instance [⁶]):

Definition 1: In input–output terms, a transcriptional module maps an excitatory input signal consisted of the concentration of a transcription factor into an output–signal consisted also of the concentration of a expressed transcription factor. If we denote z and x the genes coding the input signal and the output signal, respectively, the dynamical behavior of the transcriptional module is described in lumped time-invariant terms by the set of coupled ordinary–differential-equations given by (see Figure 2):

where:

Fig. 2 Schematic representation of a transcriptional module. The excitatory signal [z] binds to the promoter of gene x, which promotes changes in the expression of the corresponding transcription factor [x] 

Remark 3: The function f ([z]) captures the specificity of the regulatory role played by the transcription factor [z], and in general is nonlinear in nature. The regulatory function comes as a result of the interaction between [z] and the DNA–sequence which characterizes the promoter of the regulated gene x. Depending on the design, the regulatory process can involve inhibition or promotion through enzymatic dynamics.

We introduce now the class of cyclic feedback network of synthesized transcriptional modules that concerns us.

Definition 2: Cyclic feedback network of syn-thesized transcriptional modules] We define the class of n–th order cyclic feedback networks of synthesized transcriptional modules as the closed chains of transcriptional modules described as:

(3)

with:

corresponding to the i–th transcriptional module, and with:

The simple transcriptional modular structure given by (3) has been successfully applied in the construction of actual biosynthetic oscillators, with the transcription factors x_i, i = 1, . . . , n, required to display sustained oscillations through the appropriate choice of the nonlinear regulatory functions f_i, i = 1, . . . , n, and the right tuning of the parameters , , and (i.e. the nonlinear regulatory functions as well as the parameters, including the order of the system, are chosen in order to guarantee that system (3) has a stable limit cycle). The celebrated repressilator (see for instance [⁸]) is a well–known example of this family of oscillatory modular transcriptional networks governed by negative feedback, see [⁶] and [⁵].

In the sequel we describe our inverse problem, which concerns the tuning of the mathematical model consisting of a cyclic feedback modular transcription network from information provided by a synthesized system derived from the model.

3.3 Finding the Parameter Values

In the previous subsection we presented the class of modularly structured cyclic feedback systems intended to guide the conception and construction of biosynthetic oscillators based on transcriptional regulation.

Problem 1. Consider a synthesized transcriptional network to be given and resulting from a mathematical description given by system (3) (i.e. we conceive the synthesized transcriptional network to be a designed cyclic feedback system). Assume that the synthesized transcriptional net-work displays sustained oscillatory behavior. From the knowledge of the state variables {[x₁], [], [x₂], [], . . ., [x_n], []}, of the synthesized transcriptional network (i.e. the concentration of the transcription factors with their corresponding concentration of the mRNA transcripts), find a set of parameters {, , }, for i = 1, 2, . . . , n, and the parameter values of the nonlinear transcription regulatory functions f_i (·), for i = 1, 2, . . . , n, making a corresponding tuned mathematical model to closely describe the behavior of the synthesized transcriptional network.

Remark 4: We propose in the sequel a solution for Problem 1. We conceive our proposal as a proof–of–concept and because of that, and in order to simplify our exposition, we only consider here the case concerning the well-known synthetic bio-oscillator known as the repressilator.

The repressilator

In the case corresponding to the repressilator (see [8, 6]) n = 3 in (3), and the regulatory functions f₁ (·), f₂ (·) and f₃ (·), correspond to inhibition actions. In the named symmetrical case it is assumed that the transcriptional regulatory functions (all of them with negative slope) satisfy:

with the real parameter α and m defining the inhibitory regulatory functions. Moreover, in order to simplify all are assumed to be equal to 1, and it is also assumed that the decay rates and are equal to a real constant δ. Thus, the simplified symmetrical repressilator is described as follows:

(4)

(5)

(6)

Applying the Mallet–Paret and Smith theorems (see [¹⁵]), it is proved in [⁶] that this cyclic system has only one unstable equilibrium point and displays stable oscillatory behavior when the following inequality holds:

(7)

Let us now introduce our proposed tuning procedure.

4.1 Model-based On-line Parameters Estimation

Let us introduce the following assumptions (that we require in order to ensure the solvability of our systems inverse problem):

A1 The cyclic feedback system (4)–(6) is initialized to be strictly positive.

A2 The parameter vector P belongs to some neighborhood in the parameter space, such that, the system (4)–(6) exhibits sustained oscillatory behavior.

A3 The output vector Y is fully measurable.

We proceed now to expose our parameter recovering procedure.

On–line algebraic recovering of δ

This procedure consist of multiplying both sides of (11) by the time variable t, and then, integrating by parts the resulting expression. Thus, we have the expression:

(15)

By integrating once the left-hand side of the above equation we obtain the following linear expression in terms of the unknown parameter δ:

(16)

where

(17)

For sake of simplicity, we introduce the following notation:

(18)

Remark 6: Note that the system equation (16) is not well defined at time t = 0. But, for any time t after a small open time interval of the form (0, 𝜖) with 𝜖 > 0, the denominator of the equation (16) is strictly positive by assumptions A1 and A2. That is:

We focus our attention now on the computational procedure to recover the non-available regulatory signals f_i, as well as the related parameters.

Reconstructing the time variable signals f_i(y) and the unknown parameters α and m

As can be seen, the implementation of the relations in (12) needs the time derivative of the signals y_i+3 to be recovered. To overcome this issue, we introduce the following simple high–gain observer (see [⁴]). Assuming that the signals y_i are available, we can propose the following filter:

(19)

where ξ is a small positive parameter and k is a strictly positive parameter.

Proposition 1: Let y_i(t) be a scalar continuous smooth function with its corresponding time derivatives satisfying || ≤ n_i. Then, the high–gain observer recovers the time derivative of , with bounded error given by:

(20)

Then, the states converge to , when time elapses.

Proof: see Appendix A.

Consequently, the time variant regulatory signal f_i can be reconstructed by:

where are directly obtained from the filter (19). Since all the output signals, y_i+3, present oscillatory behavior (which is one of our initial assumptions), this signals and their respective time derivatives are all bounded.

Now we are able to numerically solve the one variable nonlinear algebraic equation (14) and then be able to reconstruct the parameter m. To solve this equation we can apply the standard iterative Newton-Raphson method as follows:

(21)

where z_y(m) is as defined in (14) and

Note that if , then can be computed using the formula in (13). That is:

(22)

Remark 7: Note that (21) requires to sample the outputs y_i = x_i, for i = {1, 2, 3}, since (see (14)). For this we define the indexed variables x_ik = x_i(t_k), where t_k+1 − t_k = T , with the sampling time T > 0 fixed. We compute as long as x_1k ≠ x_2k ≠ x_3k.

We proceed now to illustrate our proposal with a numerical example.

5.1 The Simulated Repressilator

We simulate the repressilator’s behavior choosing the parameter values: m = 1.5, δ = 1, and (this selection guarantees sustained oscillatory behavior). As far as the initial conditions are arbitrary, we select them as:

5.2 Algebraic Reconstruction of the Parameter δ

The numerical implementation of our reconstructor is done by using the Matlab^TM numerical computing environment. We use the Runge–Kutta method with an integration step of 0.5 × 10⁻³. The high-gain observer gains are fixed as ξ = 0.5×10⁻² and k = 1. As already established, we need to recover the parameter δ first. In Figure 3, the reconstruction of parameter δ is shown.

Fig. 3 Reconstruction error for the parameter δ 

As we can see, the parameter δ can be satis-factorily recovered regardless of the availability of measurements of , , and (in fact, notice that for the current illustrative example a set of measurements is at disposal for identification).

5.3 Estimation of the Non-available Regulatory Signals f₁, f₂ and f₃

Now, since we effectively recover δ we are able to reconstruct the non–available regulatory signals f₁, f₂ and f₃. Figure 4 shows the corresponding estimation errors for these three estimated signals. As we can see, the convergence time is very short because the high–gain observer gains were selected to be large enough.

Fig. 4 Estimation errors corresponding to the estimated signals f₁, f₂, and f₃ 

5.4 Recovering Parameters m and α

Finally, having effectively estimated the regulatory signals b f_i, the parameters m and α are effectively recovered, in that order, by solving the expression in (21) and (22) applying the Newton–Raphson method using samples of the regulatory signals b f_i and x_i every 0.5 seconds. As we can see, these parameters are in fact recovered with a very high accuracy (see Figure 5).

Fig. 5 Reconstruction parameters m and α