# Simple molecular networks that respond optimally to time-periodic stimulation

- Axel Cournac
^{1, 2}and - Jacques-Alexandre Sepulchre
^{1}Email author

**3**:29

https://doi.org/10.1186/1752-0509-3-29

© Cournac and Sepulchre; licensee BioMed Central Ltd. 2009

**Received: **11 August 2008

**Accepted: **03 March 2009

**Published: **03 March 2009

## Abstract

### Background

Bacteria or cells receive many signals from their environment and from other organisms. In order to process this large amount of information, Systems Biology shows that a central role is played by regulatory networks composed of genes and proteins. The objective of this paper is to present and to discuss simple regulatory network motifs having the property to maximize their responses under time-periodic stimulations. In elucidating the mechanisms underlying these responses through simple networks the goal is to pinpoint general principles which optimize the oscillatory responses of molecular networks.

### Results

We took a look at basic network motifs studied in the literature such as the Incoherent Feedforward Loop (IFFL) or the interlerlocked negative feedback loop. The former is also generalized to a diamond pattern, with network components being either purely genetic or combining genetic and signaling pathways. Using standard mathematics and numerical simulations, we explain the types of responses exhibited by the IFFL with respect to a train of periodic pulses. We show that this system has a non-vanishing response only if the inter-pulse interval is above a threshold. A slight generalisation of the IFFL (the diamond) is shown to work as an ideal pass-band filter. We next show a mechanism by which average of oscillatory response can be maximized by bursting temporal patterns. Finally we study the interlerlocked negative feedback loop, i.e. a 2-gene motif forming a loop where the nodes respectively activate and repress each other, and show situations where this system possesses a resonance under periodic stimulation.

### Conclusion

We present several simple motif designs of molecular networks producing optimal output in response to periodic stimulations of the system. The identified mechanisms are simple and based on known network motifs in the literature, so that that they could be embodied in existing organisms, or easily implementable by means of synthetic biology. Moreover we show that these designs can be studied in different contexts of molecular biology, as for example in genetic networks or in signaling pathways.

## Background

For past several decades, the experimental data concerning protein concentrations in biological regulatory networks have mostly been obtained in conditions related to stationary states. In some natural contexts, however, or in some experiments, the relevant cellular response to measure is the one the cell displays when it is exposed to time-dependent signals. Such information can be essential to unravelling the regulatory principles of those molecular networks which by nature are repeatedly stimulated by time-varying inputs [1]. For example the gene expression within neuronal cells can be substantially affected by the time-dependent signals received from its afferent neurons, and this property is essential for the formation of memories [2]. Other examples can be met with in various contexts, like the frequency encoding phenomenon associated with Ca^{2+} oscillations [3–6], or the response to time-dependent osmolarity shocks [7]. We note that, in these systems an optimal-response may exist with respect to periodic stimulations, meaning that for example, the mean production of some activated transcription factors, or of some proteins of interest, would be maximized if the external periodic signal follows some specific time-course and shape. This question has received little attention in the context of biological networks. This research topic is nevertheless timely, as recent developments of experimental techniques in molecular biology enables one to access time-dependent concentrations, and investigate new problems about the time-response of molecular networks.

From an experimental point of view, some interesting experiments were recently performed to reveal such time-response properties in biological systems [7–9]. In order to anticipate more experimental studies in this direction, the objective of this paper is to identify basic network topologies which allow for the property of optimal oscillatory responses in molecular biology. In the current literature only a small number of studies have tackled this issue and most of the articles addressing this question principally considers the linear response to sinusoidal perturbations. For example several papers demonstrate the low-pass filter property of small molecular networks [10–13]. Nearly all of these studies consider small amplitude perturbations, and use Fourier transform, discussing how the system parameters modify the cut-off frequency. As a matter of fact, when a system behaves like a low-pass filter, the frequency which maximizes the output response is 0. So in this case the system does not display a genuine temporal specificity. Very few studies have looked for the possibility of maximizing some variables at non-zero frequency in biological networks. A recent study [14], using again the tools of linearization and of Fourier analysis, shows that a signaling cascade with a negative feedback can behave as a band-pass filter, characterized by a frequency which maximizes the amplitude of the response. But when the system's nonlinearity is important, as is often the case in models of transcriptional regulatory networks, or when the system is linear but the periodic forcing is multiplicative (thus not additive), then Fourier decomposition of the input signal is no longer of general interest. Furthermore, in many instances, the time-dependent signals directed to regulatory networks are far from harmonic, they are often pulsatile. This is why a time-periodic input to molecular networks is often assumed as a periodic train of pulses. Smolen et al. studied a molecular network which was periodically stimulated by a train of spikes and which exhibited the ability to respond in a maximal fashion to a particular inter-spike interval [15]. Their model equations were nonlinear and associated with a lot of parameters, so the origin of the bandpass filter property is not straightforward. Another instance of periodic stimulations with a train of square pulses was thoroughly studied by Li and Goldbeter in the context of cellular receptor systems [16]. By using analytical calculations and numerical simulations of a linear model [17], they show that can exist an optimum stimulus pattern of periodic square pulses which maximizes a lumped observable. This is defined by the authors as the receptor activity. However, the Li and Goldbeter model has again many parameters, and the identification of the essential ingredients which enable one to induce a maximum amount of receptor activity is not obvious.

In the present paper we study simpler models, which pertain to basic network motifs found in biological regulations. We point out some principles which guarantee the existence of an optimal response in the production rates of output molecules, when these systems are activated by a periodic signal. The identified mechanisms are based on known network motifs in the literature and include, for example, the negative feedback loop and the incoherent feedforward loop (IFFL) motif studied by Alon and coworkers in the context of gene regulatory networks of bacteria [18]. Let us recall that in this framework a network motif is a small pattern of molecular interactions which recurs repeatedly in comparaison with what would occur in random networks. Here we will make use of this concept of motif in a broader sense, without requiring that its ubiquity be statistically proved, but requesting that it represents a basic network pattern associated with specific information-processing properties [19]. Moreover we will study these designs in at least two different contexts of molecular biology, namely in genetic regulatory networks or in signaling pathways. In the context of genetic regulation, the mathematical models are usually strongly nonlinear, involving sigmoidal regulation functions. In order to use analytical tools to study the behavior of these systems under time-dependent signals, we will employ the "logic" approximation, where regulation functions are replaced by all-or-none functions. In this case the nonlinear system is not approximated by a set of linear equations but by a set of piecewise-linear equations [20]. This approach has indeed been developed for many years in the study of genetic networks and has revealed to be quite fruitful [21].

## Results

### Simple models of biomolecule activations

In this Section, we give a brief summary of the simplest models used in the literature to describe activation processes, either at the gene level, or at the level of post-translational modification of proteins. We then report elementary but useful results which are obtained when the activation is implemented by means of a train of periodic square pulses.

*X*,

*Y*, ... The same notation is used for the concentrations of the corresponding proteins synthesized by these genes. We first introduce what is called the "simple regulation" by Alon [18], a situation which is depicted in Fig. 1(a) by an arrow from

*X*to

*Y*. In this case the gene

*Y*is positively regulated by only one transcription factor

*X*, and the basic process of synthesis of protein

*Y*from its gene is described by a single differential equation:

*R*(

*X*) is called the regulatory function. We consider regulatory functions bounded by 1. Typically

*R*(

*X*) has the form of a Hill function

*X*

^{ n }/(

*X*

^{ n }+

*θ*

^{ n }) with some cooperativity

*n*and an activation threshold

*θ*but we will often use the

*logic*approximation for

*R*(

*X*), where the latter is replaced by the unit step function

*H*(

*X*-

*θ*) (with only two values

*H*(

*X*) = 0 if X

*<*0, and

*H*(

*X*) = 1 otherwise). The parameters

*β*and

*α*are respectively the maximum synthesis rate of protein

*Y*and the degradation parameter, including the possible dilution effect from cell growth. Whenever

*R*(

*X*) = 1 the system converges to a steady state

*κ*=

*β*/

*α*. In the sequel

*X*(

*t*) is considered as a function of time, and one defines

*S*(

*t*) =

*R*(

*X*(

*t*)). A class of signals

*S*(

*t*) which will be considered below is a periodic train of square pulses of amplitude 1, whose temporal pattern is characterized by the numbers (

*τ*,

*σ*) (cf. Fig. 2). The parameter

*τ*describes the duration of the "on-phase" corresponding to the activation of transcription factor

*X*binding gene

*Y*. The inter-pulse interval, or the silent phase between the pulses, is denoted by

*σ*. Thus the period of

*S*(

*t*) is given by

*T*=

*τ*+

*σ*.

In response to the input signal *S*(*t*), we will focus on three observables in this system, namely the extrema (minimum and maximum), and the average concentrations reached by *Y* (*t*), denoted respectively by *Y*_{min}, *Y*_{max} and ⟨*Y*⟩_{
T
}. The extreme values of *Y* (*t*) can be the most relevant quantities in systems where *Y* is itself a transcription factor for other genes, because then what matters is the comparison of *Y* with a threshold of activation or of repression. On the other hand, if *Y* is consumed by some downstream process, then the mean concentration of *Y* (*t*) averaged over a period *T* of *S*(*t*), denoted by ⟨*Y*⟩_{
T
}, is the relevant quantity since it quantifies the available protein synthesized over one period.

*Y*⟩

_{ T },

*Y*

_{max}and

*Y*

_{min}as:

*σ*is varied, for a fixed pulse duration

*τ*. When

*σ*increases from 0 the observables are all decreasing. The maximum

*Y*

_{max}stays within the interval [

*Y*

_{ L },

*κ*], where

*Y*

_{ L }is the asymptotic value reached by

*Y*

_{max}(dotted line on Fig. 3, obtained when the denominator of

*Y*

_{max}, eq.(2), equals 1). We note that this level can be controlled by choosing the pulse duration

*τ*. This can be useful if

*Y*is itself a transcription factor with respect to a target gene

*Z*(see next Section).

^{2+}oscillations, whose temporal spiking patterns have been showed to significantly modulate the activity of CaM kinase II [3–6]. The activation/deactivation motif is usually represented as in Fig. 1(b) by drawing two arrows forming a cycle between two molecular states. A third arrow indicates the presence of the external signal

*S*(

*t*). Here

*Y*stands for the concentration of deactivated proteins, and

*Y** the concentration of activated ones, assuming that the total protein concentration is a conserved quantity,

*Y*+

*Y**=

*Y*

_{ tot }. Several theoretical models with different levels of complexity can be used to describe the dynamics of the interconversion between

*Y*and

*Y**. For a recent review, see e.g. [22]. The simplest form of equations which describe the dynamics of the cycle between these proteins is given in the following [23]:

where *k* and *k'* denote kinetic constants of the respectively activating and deactivating reactions. When *S*(*t*) is a periodic train of square pulses as in the above observables *Y*_{max}, *Y*_{min} and ⟨*Y*⟩_{
T
}can again be analytically calculated as a function of (*τ*, *σ*). However, their expressions are slightly more complicated than those of eqs.(2) due to the difference in time scales for the "on-phase" and "off-phase" of *Y* (*t*). For a fixed value of *τ* the graph of these functions of *σ* is qualitatively similar as seen in Fig. 3. The extrema and the mean values of *Y* (*t*) decrease when the inter-pulse interval *σ* is increased. Let us note that the same behavior is also achieved in the general case of nonlinear equations representing the dynamics of the covalent modifications by mean of the Goldbeter-Koshland model [3, 24] (cf. additional file 1).

We conclude from elementary calculations that a pulsatile activation of the simple regulation schemes (1) or (3) provides a way to adjust the variations of *Y* in a range of values which can be controlled by tuning the temporal pattern (*τ*, *σ*) of the pulses.

### Periodic activation of the IFFL

In this Section we show that one basic design which enables one to obtain an optimal output in response to a periodic train of pulses, is the *incoherent feedforward loop* (IFFL) motif studied by Alon and co-workers in the context of transcriptional networks of bacteria [25]. We first give the definition of this motif and then study new properties which appear in this system under periodic stimulations.

#### The type-I Incoherent Feedforward Loop

*X*promotes the expression of both genes

*Y*and

*Z*but the expression of the latter is repressed by

*Y*. In order to simplify the formalism so as to obtain simple analytical estimates, we again consider the logic approximation, where the regulation functions are described by step functions which take only binary values 0 or 1. As in the above,

*S*(

*t*) =

*R*(

*X*(

*t*)) denotes an input signal emanating from

*X*and the dynamics of the system is described by the following equations:

Here the first equation is a simple regulation as introduced in the previous section, and the second equation governs variable *Z* which is activated only when both the signal *S*(*t*) is "on" and when the repressor level *Y* is smaller than the threshold *θ*. We note that for simplicity's sake we have considered same degradation (*α*) and maximal production rates (*β*), for *Y* and *Z*, but this choice is not essential for the results reported below.

*Z*concentration when it is submitted to a constant stimulus. Once the signal associated to

*X*is switched "on", the variable

*Z*rises due to its positive regulation by

*X*, but this activation is terminated by the repressor

*Y*which is also activated by

*X*. The duration of the pulse can be quantified by the time

*τ*

_{ θ }needed by

*Y*to reach its repression threshold

*θ*after the onset of the signal. For the simple model above, the pulse duration is merely computed as:

This will be called the proper pulse duration of the IFFL system. Thus if the activation lasts longer than the proper pulse *τ*_{
θ
}, the system can recover its primary state. In other words the IFFL system can only detect changes of stimulation, meaning that it *adapts* to a constant stimulus. So, as mentioned by Li and Goldbeter [16] in regarding the context of receptor desensitization, in order to reach a given level of synthesis of *Z* in such system, a pulsatile pattern of stimulation must be considered instead of a continuous stimulation. Therefore in the following we consider a periodic stimulus *S*(*t*), which has the form of a square-wave, similar to that studied in the previous Section.

#### Pulsatile periodic activation of the IFFL

*Z*⟩

_{ T }and

*Z*

_{max}. As illustrated in Fig. 4 and Fig. 5, the average value ⟨

*Z*⟩

_{ T }, as well as

*Z*

_{max}can reach maximal values for specific choices of the pulse pattern (

*τ*,

*σ*). We will describe this phenomenon in more detail.

*Y*(

*t*) is identical to the one obtained in the previous section for a simple regulation, because eq.(4) coincides with eq.(1). However, in what concerns the evolution of the output protein

*Z*, a detailed analysis of this system (see additional file 1) shows that two cases should be distinguished according to the size of the pulse duration

*τ*compared with the proper pulse

*τ*

_{ θ }of the IFFL [cf. eq. (6)]. Here we detail only the case of short pulses, corresponding to the condition

*τ*<

*τ*

_{ θ }. Fig. 4 shows a 3D-plot of ⟨

*Z*⟩

_{ T }computed in function of the pulse parameters (

*τ*,

*σ*). One observes that the extrema of ⟨

*Z*⟩

_{ T }form a ridge with a roughly constant slope

*τ*/

*σ*in the parameter space (

*τ*,

*σ*). By fixing the pulse duration

*τ*Figs. 5 show the existence of maxima in ⟨

*Z*⟩

_{ T }, as well as in

*Z*

_{max}in Fig. 5b, reached for a critical inter-pulse interval

*σ*

_{1}. Another feature of these observables shown in Figs. 5 is that they start to increase from zero only if the inter-pulse interval

*σ*is above a minimum value

*σ*

_{0}. When

*τ*<

*τ*

_{ θ }, the values of

*σ*

_{0}and of

*σ*

_{1}can be analytically computed as follows (cf. Suppl. Info.):

*σ*

_{1}which should be waited for between the pulses to maximize the response of the IFFL network. We note that when

*ατ*is small, one can use the first-order approximation

*σ*

_{0}~

*σ*

_{1}~

*τ*(

*κ*/

*θ*-1), where

*κ*is the equilibrium concentration of

*Y*proteins, and

*θ*<

*κ*is the repression threshold of

*Y*with respect to gene

*Z*. Our estimate means that provided

*ατ*is small enough, the optimal pulse pattern is given by only one condition on the ratio between the on/off phases, i.e.

*σ*

_{1}/

*τ*~

*κ*/

*θ*- 1. It does not impose a unique value for the pulse duration

*τ*. Furthermore, it can be seen that the height of the optimal ⟨

*Z*⟩

_{ T }and

*Z*

_{max}is of the order of

*θ*. Thus the average and maximum productions of

*Z*are optimized when

*θ*is close to its maximal possible value, which is bounded by

*κ*. But then

*σ*

_{1}≪

*τ*, meaning that in this case the optimal rest interval would be smaller than the pulse duration. Therefore for the IFFL with a high repressor threshold, optimal trains of square pulses would seem to be alike in a constant stimulation. However, for this system, there is a crucial difference between the situation of constant stimulus, where

*Z*falls to zero, because it is completely repressed by

*Y*, and the other situation where the activation of

*Y*is periodically interrupted during for a small lapse of time. In real biochemical systems, however, the value of the threshold

*θ*cannot be arbitrarily chosen, and the implementation of our proposed principle depends on the context, being for instance the one of genetic regulation or the one of intracellular signaling. Based on plausible parameter values found in the literature, Table 1 gives examples of optimal pulse patterns predicted by our simple model in the case of genetic expression, eqs.(4)-(5), or in the case of signaling cycles, eqs.(3)–(5). The first line gives optimal interpulse durations

*σ*

_{1}for a simple IFFL composed of three genes with parameters in the range given in the experimental work of [26]. The second line gives an estimation of

*σ*

_{1}for an IFFL whose node

*Y*is a signaling cycle. The numerical values for the parameters of the signaling cycle, (

*k'*=

*k*= 1 s

^{-1},

*Y*

_{ tot }= 2

*μ*M), are based on values estimated in [27].

In conclusion one sees that the IFFL network motif (Fig. 1(c)) studied by Alon *et al*. has a nice ability to provide one of the simplest mechanisms which maximizes the average synthetized protein ⟨*Z* ⟩_{
T
}when periodically stimulated with a train of pulses.

### Mixed regulatory networks with a "Diamond" IFFL

The periodic stimulation of the IFFL network with train of pulses reveals an interesting property of optimal response with respect to the pulse pattern. In this section we slightly extend the IFFL motif to obtain a design which improves its ability of bandpass filtering, i.e. the property to induce a non zero response only for a finite range of pulse periods. We consider the "diamond IFFL" motif depicted in Fig. 1(d). The scheme is analog to the IFFL but the activation of the target gene takes two nodes instead of only one. In fact, this extension of the IFFL can be studied in the gene network seen in Fig. 1(d). It can also be implemented in a mixed regulatory network which blends signaling and gene nodes in the same interaction graph (cf. Fig. 1(e)). This graph may schematically represent a transduction network consisting of a crosstalk between two pathways converging with opposite interactions in the same genetic system. A similar scheme has been used in the study of the formation of long term memory at the molecular level [15]. The comparison between the models of the literature and the outcome of our analysis will be provided in the Discussion section.

*Z*. The equations of the system depicted on Fig. 1(e) can be written as previously with the following dynamics for the proteins ${Y}_{1}^{\ast}$, ${Y}_{2}^{\ast}$ and

*Z*:

In this simple setting where identical kinetic parameters are assumed for the covalent modification cycles 1 and 2, the time-evolutions of ${Y}_{1}^{\ast}$ and ${Y}_{2}^{\ast}$ are synchronized by the periodic signal *S*(*t*). In this case, eq. (11) shows that the expression of *Z* gene is activated only when *θ*_{1} <${Y}_{1}^{\ast}$ <*θ*_{2} (*i* = 1, 2), which entails that the relation *θ*_{1} <*θ*_{2} must be assumed in order to have a non zero production of *Z*. This assumption means that the affinity of the activator for binding the promoter of *Z* gene is higher than the affinity of the repressor for the same promoter.

*τ*,

*σ*) such that

*σ*is varying and the pulse duration

*τ*<${\tau}_{{\theta}_{1}}$ is fixed (cf. eq.(6)). Fig. 6 shows the variation of the mean value of synthesized

*Z*, denoted by ⟨

*Z*⟩

_{ T }. As for the IFFL, there is a lower limit of the inter-spike interval

*σ*

_{0}in which no production of the target gene

*Z*occurs. But once

*σ*is increased above

*σ*

_{0}there is a fast growth of ⟨

*Z*⟩

_{ T }which culminates at

*σ*

_{1}which can be computed in the same way as before (eq.(8)). If

*σ*is further augmented, ⟨

*Z*⟩

_{ T }registers a decrease which reaches 0 when

*σ*>

*σ*

_{2}on the graph. This property is due to the existence of a new node (activator Y1) in the "diamond" IFFL, which ceases to activate the target gene Z once

*Y*

_{1}max becomes smaller than its activation threshold of Z. Thus the new feature of the pulse-response of ⟨

*Z*⟩

_{ T }is that (for a fixed

*τ*) one obtains only a finite interval [

*σ*

_{0},

*σ*

_{2}] of inter-pulse intervals in which the protein

*Z*is produced. Therefore, the new property appearing in this system, as compared with the standard IFFL, gives us the possibility to use this architecture as a bandpass filtering, allowing the gene to respond only to specific periodic trains of pulses, with the inter-pulse

*σ*belonging to a limited range [

*σ*

_{0},

*σ*

_{2}].

### Optimal output in response to bursting oscillations

The network motifs studied in the previous Sections can be combined in modules that can be linked in a larger network to process complex signals. We show in this Section that by linking two diamond IFFL motifs one obtains a class of regulatory networks which are able to selectively respond to bursting oscillations.

Bursting oscillations are time-dependent signals typically emitted by neuron cells [28]. Their temporal pattern is formed by trains of spikes alternating with refractory periods. The neuronal signal is a depolarization wave created by the dynamics of ionic channels in the neuron membranes. However, in the vicinity of a synapse, this electrical signal transforms into a chemical signal affecting the receptors of the post-synaptic neuron. It thereby propagates its influence in the molecular regulatory networks of this neuron. Consequently, a bursting neuronal signal can potentiality impact on the gene expression of the post-synaptic neuron. Therefore it is worthwhile to identify simple topologies of regulatory networks with the ability to be sensitive and selective to bursting oscillations.

In order to explore this question we consider the network represented in Fig. 1(f), which combines in series two diamond IFFL motifs. Each node of Fig. 1(f) is represented by single species but in practice the proteins *Y*_{1} and *Y*_{2} on this graph could be replaced by covalent modification cycles as the ones depicted in Fig. 1(e). In this case Fig. 1(f) corresponds to a regulatory system where an input signal influences the expression of a target gene *W* by the intermediary of a signaling pathway, in which 2 phosphorylation/dephosphorylation cycles respectively activate and repress the promoter of the intermediate gene *Z*. The corresponding protein activates in turn the repressor *V*_{1} and the activator *V*_{2} of a target gene *W*.

*S*(

*t*) we consider a time-periodic pattern characterized by four time intervals (

*τ*

_{1},

*σ*

_{1},

*τ*

_{2},

*σ*

_{2}) (cf. Fig. 7). The period of the signal is given by

*T*=

*τ*

_{2}+

*σ*

_{2}, where

*τ*

_{2}is the duration of the bursting phase and

*σ*

_{2}is the quiescent period. The bursting phase is defined by a series of square spikes in which

*τ*

_{1}is the width of the spike and

*σ*

_{1}is the interspike interval. The quiescent period which separates the bursting phases is given by

*σ*

_{2}. Thus the signal is characterized by two different time-scales

*τ*

_{1}+

*σ*

_{1}≪

*τ*

_{2}+

*σ*

_{2}. In the numerical simulations reported below these time scales are separated by two orders of magnitude, reflecting the typical temporal differences between the dynamics in signaling cascades and in genetic regulations (cf. Table 1).

*S*(

*t*) described above, it has the property to selectively recognize temporal patterns of bursting oscillations, allowing the target gene to be maximally transcribed in some conditions. To illustrate this point Fig. 8 shows the time-evolution of

*W*in response to various periodic stimulations

*S*(

*t*). When the system is submitted to periodic trains of square waves without bursting (Fig. 8(a),

*σ*

_{1}= 0), the gene

*W*is not expressed. Likewise if the stimulation consists of a long train of spikes without any quiescent period (Fig. 8(b)), the average level of

*W*remains negligible. However, if we stimulate the motif with a specific bursting signal (Fig. 8(c)), the system gives a non-zero response. More generally, the striking feature of the network of Fig. 1(f) is to exhibit a non vanishing response only in a given range of pulse patterns. Moreover, if the time intervals

*τ*

_{1}and

*τ*

_{2}are fixed, the system possesses a set of maxima for some optimal values of (

*σ*

_{1},

*σ*

_{2}). In view of of Fig. 8(a–c), the system behaves as it filtered out low as well as high frequencies. But this conclusion is misleading since when high and low frequencies are mixed in the same input signal in the form of bursting oscillations, the system displays a non zero response in the evolution of

*W*, with the possibility of optimizing the average level of

*W*over one period.

The reason why the network in Fig. 1(f) manifests a large response to bursting oscillations is a direct consequence of the results discussed in the previous Sections. Let us again suppose that the active phase of the bursting lasts *τ*_{2} and the quiescent interval is *σ*_{2}. During this active phase *τ*_{2}, and assuming the spike duration *τ*_{1} is fixed, an optimal output in the level of the intermediate protein *Z* is produced if the interspike interval *σ*_{1} is chosen as the one which maximizes the response curve of Fig. 6. For example, using the parameter of that figure, this corresponds to *σ*_{1} ≃ 1/(4*α*_{1}). On the other hand, during the refractory period *σ*_{2} of the bursting oscillations, the activity of *Z* returns to zero. Thus the time-evolution of the intermediate protein *Z* is close to a periodic train of square pulses with time pattern (*τ*_{2}, *σ*_{2}) characterizing respectively its "on" and "off" phases. Now this variable feeds the regulation functions of the second subnetwork *V*_{1} - *V*_{2} - *W*. But for a fixed *τ*_{2} this subnetwork responds in an optimal way if the quiescent phase is chosen again as the critical *σ*_{2} which maximizes a response curve like the one shown in Fig. 6. For example, using the same parameter as for this figure, *σ*_{2} would correspond to 1/(4*α*_{2}). Therefore, since the signaling time scale 1/*α*_{1} is much smaller than the one of gene regulation 1/*α*_{2}, we conclude that the temporal pattern which maximizes the response in *W* is *σ*_{1} ≪ *σ*_{2}. This agrees with the concept of bursting oscillations.

### The interlocked negative feedback loop

*average*concentration of the output proteins. It is also interesting from a physiological point of view, to find biological networks that give an optimum frequency maximizing the oscillatory

*amplitude*of the output. So, we consider another basic network motif which provides a maximal response in function of a periodic stimulation. The motif network, which is shown in Fig. 1(g), is similar to the IFFL, but the repressor

*Y*in addition to being activated by X is activated by Z in a feedback loop. Furthermore, Z is not directly activated by X. This negative loop is not a gene regulatory motif in the sense defined by Alon et al., but we still call it a motif as it is the building block of several molecular biology oscillators [29]. We note that often this regulatory motif is not purely genetic, but the negative control is induced by a post-transcriptional regulation [30]. Furthermore, when the species

*Z*not only activates

*Y*but also exerts an auto-activation, the system can become a biological oscillator. In fact this positive auto-regulation is essential in order to display self-sustained oscillations. For example, minimal models of circadian oscillators have been proposed, involving the pattern of activation seen in Fig. 1(g)[31]. Likewise, minimal models for the cellular cycle have also been proposed in this way [32]. This network motif also supports a model for the oscillatory dynamics of

*p*53 [33], oscillations in Neurospora [34]. It has also been studied in a synthetic biology perspective [35]. Other examples exist (cf. Table 2), like for instance a model given by Song et al. [36] who also use the interlocked negative feedback loop to describe a model of memory formation.

Different models of the interlocked negative feedback loop. This table shows various regulators which are involved in an interlocked negative feedback loop.

*Z*(or

*Y*) is influenced by periodic variations coming from

*X*. This periodic activation can be caused by a signal associated with an autonomous oscillator, or by the outcome of a signaling pathway which is periodically stimulated, as considered in the previous Section. In all the cases we represent the periodic activation of unit

*X*in Fig. 1(g) by a signal

*S*(

*t*), and the dynamics of the protein concentrations

*Y*and

*Z*are determined by the following equations:

In this system, the signal *S*(*t*) associated with *X*(*t*) can be a periodic square pulse of amplitude 1, *R*^{+}(*Z*) is an activating regulation function, like the Hill function introduced above, and *R*^{-}(*Y*) = 1 - *R*^{+}(*Y*) is a repressing function.

The resonance response is much clearer in this case if the input signal is sinusoidal rather than a square signal. The work performed by Lipan and Wong [37] proposes the use of oscillatory signals for studying genetic networks. These authors have devised a promising experimental procedure which consists of activating and deactivating the promoter efficiency with the use of electromagnetic fields. The net effect of their procedure amounts to periodically modulating the expression rate of genes.

*β*

_{1}= 0, the system possesses a non-zero steady state (

*Y*

_{ s },

*Z*

_{ s }) which can be shown to always be a stable focus. The corresponding frequency

*ω*

_{0}is easily determined by computing the imaginary part of the Jacobian eigenvalues. Thus, if the system is drawn away from equilibrium, it will return to its steady state by exhibiting transient oscillations with the frequency

*ω*

_{0}. This frequency can also be revealed by periodically modulating the expression rate of

*Y*. In this case we suppose that the modulation is operated with a frequency varying in a range around

*ω*

_{0}. Numerical simulations of eqs. (18) show that in this range the network dynamics give rise to a resonance which coincides with

*ω*

_{0}. In Fig. 9, we can see a maximum in the amplitude response of protein Y. A typical order or magnitude for the resonant period obtained from plausible parameters is

*T*

_{0}~300 min. Moreover we observe that the resonance width mainly depends on the stiffness of the regulatory functions, quantified here by the cooperativity coefficient

*n*of the Hill function. To explain this property, an easy but instructive example which can be analytically dealt with, considers identical thresholds

*θ*for the regulatory Hill functions

*R*

^{+}(

*Z*) and

*R*

^{-}(

*Y*), and also assume

*β*

_{0}=

*β*. The steady state for

*Y*and

*Z*is symmetric, i.e.

*Y*

_{ s }=

*Z*

_{ s }=

*κ*/2 (we recall that

*κ*=

*β*/

*α*). The natural frequency can be written as

*ω*

_{0}=

*nα*/2, where

*n*is the Hill coefficient of functions

*R*

^{+}(

*Z*) and

*R*

^{-}(

*Y*). Since, in absence of periodic stimulations, the damping of the oscillations is characterized by

*α*, we can estimate the "quality" factor of this resonance with the ratio

*ω*

_{0}/

*α*=

*n*/2. This indeed depends only on the Hill coefficient

*n*. Examples of resonances with

*n*= 2 and

*n*= 4 are illustrated in Fig. 9.

## Discussion

In this study we have looked at simple molecular networks showing the ability to maximize their responses to periodic stimulations. Using analytical and numerical tools, we have identified two types of basic network motifs, and some of their straightforward extensions, which possess this property (Fig. 1). The first of thes motifs is the Incoherent Feedforward Loop network (IFFL), for which we have shown that periodic activations with a train of square pulses lead to maximizing the time-average production of the output protein over one period provided that the pulse and the inter-pulse intervals are adequately chosen (Fig. 4). When the pulse duration is short enough, this optimization also applies to the maximal concentration of the output protein reached over one period (Fig. 5b). The mathematical analysis of this optimization property is simple, although it does not reckon as do most current studies which focus on this topic, with the assumption of small amplitude stimulations, which allows to work in the framework of linear filter analysis. We have also pointed out another network motif, namely the interlocked negative feedback loop which has the property to exhibit a maximum in the amplitude of oscillations when it is submitted to periodic stimulations. In this case the maximization is observed in the linear regime (i.e. for small amplitude of the stimulations) and it can be described as a standard phenomenon of resonance.

The interlocked negative loop is a well-studied motif in the context of oscillating regulatory networks (cf. examples given on Table 2). In all these examples the auto-activation of one of the two nodes of the network is imperative in order to produce sustained oscillations. Here we pinpoint a property of the interlocked feedback loop that exists without or with only a weak auto-activation. In this case, the system relaxes towards a stationary state in absence of periodic stimuli. It nevertheless possesses the ability to show up a resonant behavior which can be waked by periodic stimulations. This phenomenon might be important in some situations. For example the protein P53 is engaged in many regulation processes and it is known to oscillate with a period of 5.5 hours in some conditions. In [38] the authors show that the amplitude of oscillations of P53 is variable whereas the period of oscillations has a spectacular regularity. We speculate that even in absence of autonomous oscillations this system might show up an amplified response to external periodic stimulations (e.g. a periodic exposure to UV rays) because of the resonance phenomenon. A similar remark can be made regarding regulatory systems related to stress responses as in the SOS regulator LexA of *E. coli* or in the NF-*κ* B system.

As an illustrative application of the periodic activation of the diamond IFFL motif we believe that the feature studied in the present paper reveals the principle behind an intriguing phenomenon studied by Smolen et al. [15]. In this reference, the authors proposed a model for the formation of long term memory (LTM), a process which is based on the LTP phenomenon (Long Term Potentiation) in neurons and on the subsequent strengthening of synapses. The LTP appears in a synapse when the postsynaptic cell is submitted to some external signal, for instance in the form of glutamate pulses of significant amplitude. Then this signal induces the phosphorylation of transcription factors which genetically activate the production of new proteins in order to strengthen connections between synapses or to create new synapses. On the other hand, it has been shown that in order to learn, several animals or organisms need to receive stimulations in repeated sessions (a well-studied example is the aplysia, see [39] for a review). Interestingly these sessions must be spaced with a minimal interval. Massed training produced significantly shorter lasting memory as experimentally shown in [40]. Thus the genetic network underlying this regulation must respond only if the learning sessions are sufficiently spaced.

The model of Smolen et al. (Fig. 5 of [15]) had the same structure as the diamond IFFL shown in Fig. 1(e). It postulated the existence of two transcriptional factors, respectively an activator and a repressor, whose phosphorylated forms are activated by an external signal, allowing them to act on the promoter of a gene involved in the formation of LTM (e.g. the *creb* genes). The purpose of this modeling was to show the marked differences in the transcription rate of the LTM genes between the case of "massed stimuli" and the case of "spaced stimuli", as was experimentally observed e.g. by Yin et al. [2]. In fact an outstanding property of their model was the existence of a maximum in the transcription rate of the LTM genes, in function of the interspike interval (here corresponding to the rest time between the learning sessions). However the definite mechanism underlying this property is not easy to identify in that system, as the model counted a half dozen of nonlinear equations, with nearly 20 parameters to be chosen. Despite this complexity the authors did conclude by suggesting that a general mechanism enabling the tuning of the response of a gene promoter to an optimal stimulus frequency could be the presence of two competing processes with different sensitivities acting on this gene. In this article we have achieved a minimal model which puts this intuition on firm grounds. Our analysis shows that the presence of two competing processes related by the IFFL mechanism is indeed a key ingredient, a higher affinity of the activator compared with the one of the repressor for the target gene promoter is another one.

The results of Li and Goldbeter [16] mentioned in the introduction are another rare example found in the literature where a maximization process of the response in a molecular network is shown to critically depend on some temporal patterns of activation. Although this study is based on analytical calculations the basic ingredient which explains the observed property is not transparent. Although the link of this model with the IFFL motif is not established we are currently working on a generalization of the optimization principle identified in the IFFL responses in order to place the two studies in the same framework.

Based on typical parameter values met in different systems (genetic, signaling), Table 1 provides some estimates of the pulse and of the inter-pulse intervals which allow for maximization of the average production of target proteins. Since these estimates seem compatible with plausible experimental values, they prompt experimental implementations of the scheme. For example, several IFFL regulatory motifs have been recognized in a genome map of *E. coli* [18]). It would be interesting to stimulate one of these loops with trains of periodic pulses in order to discover whether the predicted maximization is feasible in practice. Utilization of thermosensitive promotors may be one possibility for the implementation of the periodical stimulations [41]. Nevertheless we note that the phenomena presented here may be sensitive to molecular noise and to environmental perturbations. Presumably, the envisaged experiments might be less difficult to achieve in eukaryote organisms where the genetic regulations involve more proteins than in prokaryotes, in this way, molecular noise could be reduced.

A general issue concerning the concept of network motifs, like the IFFL, is that such a module is in fact embedded in larger and more complex regulatory networks [42]. This weakens the actual meaning of modular subnetworks. In this regard, however, by considering time dependent signals we propose that frequency selectivity mechanisms like the ones studied here, can give more sense to the modular approach. Finally we note that from a biomedical perspective the question of optimal response in molecular networks could be especially relevant in situations where a maximal protein synthesis would be beneficial in periodic stimulation-based therapies. For example some innovative treatments of neural based diseases like Parkinson's, or of severe depressions, use a novel technique consisting of deep brain stimulations [43]. The patient wears a pacemaker linked to electrodes implanted in some brain nuclei. Electrical stimulations are then transmitted at a given frequency. Although this type of clinical treatment has been reported to lead to spectacular relief, the question of fine tuning the frequency of the pacemaker to achieve an optimal result is empirically solved for the time being [44]. A better understanding of the time-dependent responses of intracellular signaling pathways could be very relevant in these situations. In a different context, the traditional periodic intake of medicines is another example where optimization will emerge from rational thinking in terms of optimal pulsatile stimulations.

## Conclusion

The mechanisms presented in this article identify new possible strategies which can be employed by a cell to integrate time-dependent information provided by the environment. The studied network motifs offer the attractive possibility of selecting a signal according to its temporal structure or to its frequency content. Only messages which for instance, which have a minimal inter-pulse duration are allowed to pass. The identified mechanisms are simple and based on known network motifs in the literature, they can be embodied in existing organisms. They could also be implemented by means of synthetic biology [41, 45].

## Methods

### Numerical simulations

Numerical simulations were made with programs written in C language and the Runge-Kutta routines from the GNU Scientific Library. The simulation codes of the models which were used to obtain the numerical results reported in the paper (e.g. Figs. 5, 6, 8 and 9) are provided in the additional file 2 (with the name "archive-codes").

## Declarations

### Acknowledgements

A.C. benefits from a PhD Grant from the French Minister of Research.

## Authors’ Affiliations

## References

- Goldbeter A: Biochemical Oscillations and Cellular Rhythms: The molecular bases of periodic and chaotic behaviour. 1996, Cambridge University PressView ArticleGoogle Scholar
- Yin JC, Vecchio MD, Zhou H, Tully T: CREB as a memory modulator: induced expression of a dCREB2 activator isoform enhances long-term memory in Drosophila. Cell. 1995, 81: 107-115. 10.1016/0092-8674(95)90375-5View ArticlePubMedGoogle Scholar
- Goldbeter A, Dupont G, Berridge MJ: Minimal model for signal-induced Ca2+ oscillations and for their frequency encoding through protein phosphorylation. Proc Natl Acad Sci USA. 1990, 87 (4): 1461-1465. 10.1073/pnas.87.4.1461PubMed CentralView ArticlePubMedGoogle Scholar
- Dupont G, Goldbeter A: Protein phosphorylation driven by intracellular calcium oscillations: a kinetic analysis. Biophys Chem. 1992, 42 (3): 257-270. 10.1016/0301-4622(92)80018-ZView ArticlePubMedGoogle Scholar
- Koninck PD, Schulman H: Sensitivity of CaM kinase II to the frequency of Ca2+ oscillations. Science. 1998, 279 (5348): 227-230. 10.1126/science.279.5348.227View ArticlePubMedGoogle Scholar
- Dupont G, Houart G, Koninck PD: Sensitivity of CaM kinase II to the frequency of Ca2+ oscillations: a simple model. Cell Calcium. 2003, 34 (6): 485-497. 10.1016/S0143-4160(03)00152-0View ArticlePubMedGoogle Scholar
- Mettetal JT, Muzzey D, Gómez-Uribe C, van Oudenaarden A: The frequency dependence of osmo-adaptation in Saccharomyces cerevisiae. Science. 2008, 319 (5862): 482-484. 10.1126/science.1151582PubMed CentralView ArticlePubMedGoogle Scholar
- Dolmetsch RE, Xu K, Lewis RS: Calcium oscillations increase the efficiency and specificity of gene expression. Nature. 1998, 392 (6679): 933-936. 10.1038/31960View ArticlePubMedGoogle Scholar
- Hersen P, McClean MN, Mahadevan L, Ramanathan S: Signal processing by the HOG MAP kinase pathway. Proc Natl Acad Sci USA. 2008, 105 (20): 7165-7170. 10.1073/pnas.0710770105PubMed CentralView ArticlePubMedGoogle Scholar
- Gomez-Uribe C, Verghese GC, Mirny LA: Operating Regimes of Signaling Cycles: Statics, Dynamics, and Noise Filtering. PLoS Comput Biol. 2007, 3 (12): e246- 10.1371/journal.pcbi.0030246PubMed CentralView ArticlePubMedGoogle Scholar
- Cox CD, McCollum JM, Austin DW, Allen MS, Dar RD, Simpson ML: Frequency domain analysis of noise in simple gene circuits. Chaos. 2006, 16 (2): 026102- 10.1063/1.2204354View ArticlePubMedGoogle Scholar
- Shankaran H, Wiley HS, Resat H: Receptor downregulation and desensitization enhance the information processing ability of signalling receptors. BMC Syst Biol. 2007, 1: 48- 10.1186/1752-0509-1-48PubMed CentralView ArticlePubMedGoogle Scholar
- Shankaran H, Resat H, Wiley HS: Cell surface receptors for signal transduction and ligand transport: a design principles study. PLoS Comput Biol. 2007, 3 (6): e101- 10.1371/journal.pcbi.0030101PubMed CentralView ArticlePubMedGoogle Scholar
- Locasale JW: Signal duration and the time scale dependence of signal integration in biochemical pathways. BMC Syst Biol. 2008, 2 (1): 108-ArXiv:0802.2683v1 [q-bio.MN] 10.1186/1752-0509-2-108PubMed CentralView ArticlePubMedGoogle Scholar
- Smolen P, Baxter DA, Byrne JH: Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory systems. Am J Physiol. 1998, 274 (2 Pt 1): C531-C542.PubMedGoogle Scholar
- Li Y, Goldbeter A: Frequency specificity in intercellular communication. Influence of patterns of periodic signaling on target cell responsiveness. Biophys J. 1989, 55: 125-145. 10.1016/S0006-3495(89)82785-7PubMed CentralView ArticlePubMedGoogle Scholar
- Segel LA, Goldbeter A, Devreotes PN, Knox BE: A mechanism for exact sensory adaptation based on receptor modification. J Theor Biol. 1986, 120 (2): 151-179. 10.1016/S0022-5193(86)80171-0View ArticlePubMedGoogle Scholar
- Alon U: An Introduction to Systems Biology. 2006, Chapman & Hall/CRCGoogle Scholar
- Wolf DM, Arkin AP: Motifs, modules and games in bacteria. Curr Opin Microbiol. 2003, 6 (2): 125-134. 10.1016/S1369-5274(03)00033-XView ArticlePubMedGoogle Scholar
- Gouzé J, Sari T: A class of piecewise linear differential equations arising in biological models. Dyn Syst. 2002, 17: 299-316.View ArticleGoogle Scholar
- Glass L, Kauffman S: The logical analysis of continuous, non-linear biochemical control networks. J Theor Bio. 1973, 39: 103-129. 10.1016/0022-5193(73)90208-7.View ArticleGoogle Scholar
- Ventura AC, Sepulchre JA, Merajver SD: A hidden feedback in signaling cascades is revealed. PLoS Comput Biol. 2008, 4 (3): e1000041- 10.1371/journal.pcbi.1000041PubMed CentralView ArticlePubMedGoogle Scholar
- Heinrich R, Neel BG, Rapoport TA: Mathematical models of protein kinase signal transduction. Mol Cell. 2002, 9 (5): 957-970. 10.1016/S1097-2765(02)00528-2View ArticlePubMedGoogle Scholar
- Goldbeter A, Koshland DE: An amplified sensitivity arising from covalent modification in biological systems. Proc Natl Acad Sci USA. 1981, 78 (11): 6840-6844. 10.1073/pnas.78.11.6840PubMed CentralView ArticlePubMedGoogle Scholar
- Mangan S, Itzkovitz S, Zaslaver A, Alon U: The incoherent feed-forward loop accelerates the response-time of the gal system of Escherichia coli. J Mol Biol. 2006, 356 (5): 1073-1081. 10.1016/j.jmb.2005.12.003View ArticlePubMedGoogle Scholar
- Rosenfeld N, Young JW, Alon U, Swain PS, Elowitz MB: Gene regulation at the single-cell level. Science. 2005, 307 (5717): 1962-1965. 10.1126/science.1106914View ArticlePubMedGoogle Scholar
- Kholodenko BN: Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem. 2000, 267 (6): 1583-1588. 10.1046/j.1432-1327.2000.01197.xView ArticlePubMedGoogle Scholar
- Izhikevich EM: Bursting. Scholarpedia. 2006, 1 (3): 1300-View ArticleGoogle Scholar
- Barkai N, Leibler S: Circadian clocks limited by noise. Nature. 2000, 403 (6767): 267-268.PubMedGoogle Scholar
- Guantes R, Poyatos JF: Dynamical principles of two-component genetic oscillators. PLoS Comput Biol. 2006, 2 (3): e30- 10.1371/journal.pcbi.0020030PubMed CentralView ArticlePubMedGoogle Scholar
- Gonze D, Halloy J, Goldbeter A: Robustness of circadian rhythms with respect to molecular noise. Proc Natl Acad Sci USA. 2002, 99 (2): 673-678. 10.1073/pnas.022628299PubMed CentralView ArticlePubMedGoogle Scholar
- Tyson JJ, Csikasz-Nagy A, Novak B: The dynamics of cell cycle regulation. Bioessays. 2002, 24 (12): 1095-1109. 10.1002/bies.10191View ArticlePubMedGoogle Scholar
- Ciliberto A, Novak B, Tyson JJ: Steady states and oscillations in the p53/Mdm2 network. Cell Cycle. 2005, 4 (3): 488-493.View ArticlePubMedGoogle Scholar
- Dunlap JC, Loros JJ: The neurospora circadian system. J Biol Rhythms. 2004, 19 (5): 414-424. 10.1177/0748730404269116View ArticlePubMedGoogle Scholar
- Atkinson MR, Savageau MA, Myers JT, Ninfa AJ: Development of genetic circuitry exhibiting toggle switch or oscillatory behavior in Escherichia coli. Cell. 2003, 113 (5): 597-607. 10.1016/S0092-8674(03)00346-5View ArticlePubMedGoogle Scholar
- Song H, Smolen P, Av-Ron E, Baxter DA, Byrne JH: Dynamics of a minimal model of interlocked positive and negative feedback loops of transcriptional regulation by cAMP-response element binding proteins. Biophys J. 2007, 92 (10): 3407-3424. 10.1529/biophysj.106.096891PubMed CentralView ArticlePubMedGoogle Scholar
- Lipan O, Wong WH: The use of oscillatory signals in the study of genetic networks. Proc Natl Acad Sci USA. 2005, 102 (20): 7063-7068. 10.1073/pnas.0403790102PubMed CentralView ArticlePubMedGoogle Scholar
- Geva-Zatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, Dekel E, Yarnitzky T, Liron Y, Polak P, Lahav G, Alon U: Oscillations and variability in the p53 system. Mol Syst Biol. 2006, 2: 2006.0033- 10.1038/msb4100068PubMed CentralView ArticlePubMedGoogle Scholar
- Kandel ER: The molecular biology of memory storage: a dialog between genes and synapses. Biosci Rep. 2004, 24 (4–5): 475-522. 10.1007/s10540-005-2742-7View ArticlePubMedGoogle Scholar
- Carew TJ, Pinsker HM, Kandel ER: Long-Term Habituation of a Defensive Withdrawal Reflex in Aplysia. Science. 1972, 175 (4020): 451-454. 10.1126/science.175.4020.451View ArticlePubMedGoogle Scholar
- Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339-342. 10.1038/35002131View ArticlePubMedGoogle Scholar
- Mazurie A, Bottani S, Vergassola M: An evolutionary and functional assessment of regulatory network motifs. Genome Biol. 2005, 6 (4): R35- 10.1186/gb-2005-6-4-r35PubMed CentralView ArticlePubMedGoogle Scholar
- Perlmutter JS, Mink JW: Deep brain stimulation. Annu Rev Neurosci. 2006, 29: 229-257. 10.1146/annurev.neuro.29.051605.112824PubMed CentralView ArticlePubMedGoogle Scholar
- Henning J, Koczan D, Glass A, Karopka T, Pahnke J, Rolfs A, Benecke R, Gimsa U: Deep brain stimulation in a rat model modulates TH, CaMKIIa and Homer1 gene expression. Eur J Neurosci. 2007, 25: 239-250. 10.1111/j.1460-9568.2006.05264.xView ArticlePubMedGoogle Scholar
- Elowitz M, Leibler S: A synthetic oscillatory network of transcriptionnal regulators. Nature. 2000, 403: 335-338. 10.1038/35002125View ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.