Global entrainment of transcriptional systems to periodic inputs

This paper addresses the problem of giving conditions for transcriptional systems to be globally entrained to external periodic inputs. By using contraction theory, a powerful tool from dynamical systems theory, it is shown that certain systems driven by external periodic signals have the property that all solutions converge to a fixed limit cycle. General results are proved, and the properties are verified in the specific case of some models of transcriptional systems. The basic mathematical results needed from contraction theory are proved in the paper, making it self-contained.


Introduction
Periodic, clock-like rhythms pervade nature and regulate the function of all living organisms. For instance, circadian rhythms are regulated by an endogenous biological clock entrained by the light signals from the environment that then acts as a pacemaker [1]. Moreover, such an entrainment can be obtained even if daily variations are present, like e.g. temperature and light variations. Another important example of entrainment in biological systems is at the molecular level, where the synchronization of several cellular processes is regulated by the cell cycle [2].
An important question in mathematical and computational biology is that of finding conditions ensuring that entrainment occurs. The objective is to identify classes of biological systems that can be entrained by an exogenous signal. To solve this problem, modelers often resort to simulations in order to show the existence of periodic solutions in the system of interest. Simulations, however, can never prove that solutions will exist for all parameter values, and they are subject to numerical errors. Moreover, robustness of entrained solutions needs to be checked in the presence of noise and uncertainties, which cannot be avoided experimentally.
From a mathematical viewpoint, the problem of formally showing that entrainment takes place is known to be extremely difficult. Indeed, if a stable linear time-invariant model is used to represent the system of interest, then entrainment is usually expected, when the system is driven by an external periodic input, with the system response being a filtered, shifted version of the external driving signal. However, in general, as is often the case in biology, models are nonlinear. The response of nonlinear systems to periodic inputs is the subject of much current systems biology experimentation; for example, in [3], the case of a cell signaling system driven by a periodic square-wave input is considered. From measurements of a periodic output, the authors fit a transfer function to the system, implicitly modeling the system as linear even though (as stated in the Supplemental Materials to [3]) there are saturation effects so the true system is nonlinear. For nonlinear systems, driving the system by an external periodic signal does not guarantee the system response to also be a periodic solution, as nonlinear systems can exhibit harmonic generation or suppression and complex behavior such as chaos or quasi-periodic solutions [4]. This may happen even if the system is well-behaved with respect to constant inputs; for example, there are systems which converge to a fixed steady state no matter what is the input excitation, so long as this input signal is constant, yet respond chaotically to the simplest oscillatory input; we outline such an example in the Materials and Methods Section, see also [5]. Thus, a most interesting open problem is that of finding conditions for the entrainment to external inputs of biological systems modeled by sets of nonlinear differential equations.
One approach to analyzing the convergence behavior of nonlinear dynamical systems is to use Lyapunov functions. However, in biological applications, the appropriate Lyapunov functions are not always easy to find and, moreover, convergence is not guaranteed in general in the presence of noise and/or uncertainties. Also, such an approach can be hard to apply to the case of non-autonomous systems (that is, dynamical systems directly dependent on time), as is the case when dealing with periodically forced systems.
The above limitations can be overcome if the convergence problem is interpreted as a property of all trajectories, asking that all solutions converge towards one another (contraction). This is the viewpoint of contraction theory, [6], [7], and more generally incremental stability methods [8]. Global results are possible, and these are robust to noise, in the sense that, if a system satisfies a contraction property then trajectories remain bounded in the phase space [9]. Contraction theory has a long history. Contractions in metric functional spaces can be traced back to the work of Banach and Caccioppoli [10] and, in the field of dynamical systems, to [11] and even to [12] (see also [13], [8], and e.g. [14] for a more exhaustive list of related references). Contraction theory has been successfully applied to both nonlinear control and observer problems, [7], [15] and, more recently, to synchronization and consensus problems in complex networks [16], [17], [18]. In [19] it was proposed that contraction can be particularly useful when dealing with the analysis and characterization of biological networks. In particular, it was found that using non Euclidean norms, as also suggested in [6] (Sec. 3.7ii), can be particularly effective in this context [19], [20].
One of the objectives of this paper is to give a self-contained exposition, with all proofs included, of results in contraction theory as applied to entrainment of periodic signals, and, moreover, to show their applicability to problems of biological interest. We believe that contraction analysis should be recognized as an important component of the ''toolkit'' of systems biology, and this paper should be useful to other researchers contemplating the use of these tools.
For concreteness, we focus mainly on transcriptional systems, as well as related biochemical systems, which are basic building blocks for more complex biochemical systems. However, the results that we obtain are of more generality. To illustrate this generality, and to emphasize the use of our techniques in synthetic biology design, we discuss as well the entrainment of a Repressilator circuit in a parameter regime in which endogenous oscillations do not occur, as well as the synchronization of a network of Repressilators. A surprising fact is that, for these applications, and contrary to many engineering applications, norms other than Euclidean, and associated matrix measures, must be considered.

Mathematical tools
We consider in this paper systems of ordinary differential equations, generally time-dependent: defined for t [ ½0,?) and x [ C, where C is a subset of R n . It will be assumed that f (t,x) is differentiable on x, and that f (t,x), as well as the Jacobian of f with respect to x, denoted as both continuous in (t,x). In applications of the theory, it is often the case that C will be a closed set, for example given by non-negativity constraints on variables as well as linear equalities representing mass-conservation laws. For a non-open set C, differentiability in x means that the vector field f (t,.) can be extended as a differentiable function to some open set which includes C, and the continuity hypotheses with respect to (t,x) hold on this open set. We denote by Q(t,s,j) the value of the solution x(t) at time t of the differential equation (1) with initial value x(s)~j. It is implicit in the notation that Q(t,s,j) [ C (''forward invariance'' of the state set C). This solution is in principle defined only on some interval sƒtvsze, but we will assume that Q(t,s,j) is defined for all t §s.
Conditions which guarantee such a ''forward-completeness'' property are often satisfied in biological applications, for example whenever the set C is closed and bounded, or whenever the vector field f is bounded. (See Appendix C in [21] for more discussion, as well as [22] for a characterization of the forward completeness property.) Under the stated assumptions, the function Q is jointly differentiable in all its arguments (this is a standard fact on well-posedness of differential equations, see for example Appendix C in [21]).
We recall (see for instance [23]) that, given a vector norm on Euclidean space ( . j j), with its induced matrix norm A k k, the associated matrix measure m is defined as the directional derivative of the matrix norm, that is, For example, if . j j is the standard Euclidean 2-norm, then m(A) is the maximum eigenvalue of the symmetric part of A. As we shall see, however, different norms will be useful for our applications. Matrix measures are also known as ''logarithmic norms'', a concept independently introduced by Germund Dahlquist and Sergei Lozinskii in 1959, [24,25]. The limit is known to exist, and the convergence is monotonic, see [24,26].
We will say that system (1) is infinitesimally contracting on a convex set C(R n if there exists some norm in C, with associated matrix measure m such that, for some constant c [ R{ 0 f g, Let us discuss informally (rigorous proofs are given later) the motivation for this concept. Since by assumption f t,x ð Þ is

Author Summary
The activities of living organisms are governed by complex sets of biochemical reactions. Often, entrainment to certain external signals helps control the timing and sequencing of reactions. An important open problem is to understand the onset of entrainment and under what conditions it can be ensured in the presence of uncertainties, noise, and environmental variations. In this paper, we focus mainly on transcriptional systems, modeled by Ordinary Differential Equations. These are basic building blocks for more complex biochemical systems. However, the results that we obtain are of more generality. To illustrate this generality, and to emphasize the use of our techniques in synthetic biology, we discuss the entrainment of a Repressilator circuit and the synchronization of a network of Repressilators. We answer the following two questions: 1) What are the dynamical mechanisms that ensure the entrainment to periodic inputs in transcriptional modules? 2) Starting from natural systems, what properties can be used to design novel synthetic biological circuits that can be entrained? For some biological systems which are always ''in contact'' with a continuously changing environment, entrainment may be a ''desired'' property. Thus, answering the above two questions is of fundamental importance. While entrainment may appear obvious at first thought, it is not a generic property of nonlinear dynamical systems. The main result of our paper shows that, even if the transcriptional modules are modeled by nonlinear ODEs, they can be entrained by any (positive) periodic signal. Surprisingly, such a property is preserved if the system parameters are varied: entrainment is obtained independently of the particular biochemical conditions. We prove that combinations of the above transcriptional module also show the same property. Finally, we show how the developed tools can be applied to design synthetic biochemical systems guaranteed to exhibit entrainment.
continuously differentiable, the following exact differential relation can be obtained from (1): where, as before, J~J t,x ð Þ denotes the Jacobian of the vector field f , as a function of x [ C and t [ R z , and where Dx denotes a small change in states and ''D _ x x'' means dDx=dt, evaluated along a trajectory. (In mechanics, as in [27], Dx is called ''virtual displacement'', and formally it may be thought of as a linear tangent differential form, differentiable with respect to time.) Consider now two neighboring trajectories of (1), evolving in C, and the virtual displacements between them. Note that (3) can be thought of as a linear time-varying dynamical system of the form: once that J(t)~J(t,x(t)) is thought of as a fixed function of time.
Hence, an upper bound for the magnitude of its solutions can be obtained by means of the Coppel inequality [28], yielding: where m J ð Þ is the matrix measure of the system Jacobian induced by the norm being considered on the states and Dx 0 ð Þ j j~Dx 0 j j. Using (4) and (2), we have that Thus, trajectories starting from infinitesimally close initial conditions converge exponentially towards each other. In what follows we will refer to c 2 as contraction (or convergence) rate.
The key theoretical result about contracting systems links infinitesimal and global contractivity, and is stated below. This result can be traced, under different technical assumptions, to e.g. [6], [13], [12], [11]. Theorem 1. Suppose that C is a convex subset of R n and that f (t,x) is infinitesimally contracting with contraction rate c 2 . Then, for every two solutions x(t)~Q(t,0,j) and z(t)~Q(t,0,f) of (1), it holds that: In other words, infinitesimal contractivity implies global contractivity. In the Materials and Methods section, we provide a self-contained proof of Theorem 1. In fact, the result is shown there in a generalized form, in which convexity is replaced by a weaker constraint on the geometry of the space.
In actual applications, often one is given a system which depends implicitly on the time, t, by means of a continuous function u t ð Þ, i.e. systems dynamics are represented by _ x x~f x,u t ð Þ ð Þ. In this case, u t ð Þ : R z ?U (where U is some subset of R), represents an external input. It is important to observe that the contractivity property does not require any prior information about this external input. In fact, since u t ð Þ does not depend on the system state variables, when checking the property, it may be viewed as a constant parameter, u [ U. Thus, if contractivity of f x,u ð Þ holds uniformly Vu [ U, then it will also hold for f x,u t ð Þ ð Þ. Given a number Tw0, we will say that system (1) is T-periodic if it holds that Notice that the system _ x x~f x,u t ð Þ ð Þis T-periodic, if the external input, u t ð Þ, is itself a periodic function of period T. The following is the basic theoretical result about periodic orbits that will be used in the paper. A proof may be found in [6], Sec. 3.7.vi.
Theorem 2. Suppose that: N C is a closed convex subset of R n ; N f is infinitesimally contracting with contraction rate c 2 ; N f is T-periodic.
Then, there is a unique periodic solution a(t) : ½0,?)?C of (1) of period T and, for every solution x(t), it holds that x t ð Þ{a t ð Þ j j ?0 as t??.
In the Materials and Methods section of this paper, we provide a self-contained proof of Theorem 2, in a generalized form which does not require convexity.

A simple example
As a first example to illustrate the application of the concepts introduced so far, we choose a simple bimolecular reaction, in which a molecule of A and one of B can reversibly combine to produce a molecule of C.
This system can be modeled by the following set of differential equations: where we are using A~A(t) to denote the concentration of A and so forth. The system evolves in the positive orthant of R 3 . Solutions satisfy (stoichiometry) constraints: for some constants a and b.
We will assume that one or both of the ''kinetic constants'' k i are time-varying, with period T. Such a situation arises when the k i 's depend on concentrations of additional enzymes, which are available in large amounts compared to the concentrations of A,B,C, but whose concentrations are periodically varying. The only assumption will be that k 1 (t) §k 0 Because of the conservation laws (7), we may restrict our study to the equation for C. Once that all solutions of this equation are shown to globally converge to a periodic orbit, the same will follow for A(t)~a{C(t) and B(t)~b{C(t). We have that: Because A(t) §0 and B(t) §0, this system is studied on the subset of R defined by 0ƒCƒ min a,b f g. The equation can be rewritten as: Differentiation with respect to C of the right-hand side in the above system yields this (1|1) Jacobian: Since we know that {azCƒ0 and {bzCƒ0, it follows that Using any norm (this example is in dimension one) we have that m(J)v{c 2 . So (6) is contracting and, by means of Theorem 2, solutions will globally converge to a unique solution of period T (notice that such a solution depends on system parameters). Figure 1 shows the behavior of the dynamical system (9), using two different values of k {1 . Notice that the asymptotic behavior of the system depends on the particular choice of the biochemical parameters being used. Furthermore, it is worth noticing here that the higher the value of k {1 , the faster will be the convergence to the attractor.

Mathematical model and problem statement
We study a general externally-driven transcriptional module. We assume that the rate of production of a transcription factor X is proportional to the value of a time dependent input function u(t), and X is subject to degradation and/or dilution at a linear rate. (Later, we generalize the model to also allow nonlinear degradation as well.) The signal u(t) might be an external input, or it might represent the concentration of an enzyme or of a second messenger that activates X . In turn, X drives a downstream transcriptional module by binding to a promoter (or substrate), denoted by e with concentration e~e(t). The binding reaction of X with e is reversible and given by: where Y is the complex protein-promoter, and the binding and dissociation rates are k 1 and k 2 respectively. As the promoter is not subject to decay, its total concentration, e T , is conserved, so that the following conservation relation holds: We wish to study the behavior of solutions of the system that couples X and e, and specifically to show that, when the input u(t) is periodic with period T, this coupled system has the property that all solutions converge to some globally attracting limit cycle whose period is also T. Such transcriptional modules are ubiquitous in biology, natural as well as synthetic, and their behavior was recently studied in [29]  in the context of ''retroactivity'' (impedance or load) effects. If we think of u(t) as the concentration of a protein Z that is a transcription factor for X , and we ignore fast mRNA dynamics, such a system can be schematically represented as in Figure 2, which is adapted from [29]. Notice that u(t) here does not need to be the concentration of a transcriptional activator of X for our results to hold. The results will be valid for any mathematical model for the concentrations, x, of X and y, of Y (the concentration of e is conserved) of the form: An objective in this paper is, thus, to show that, when u is a periodic input, all solutions of system (12) converge to a (unique) limit cycle ( Figure 3). The key tool in this analysis is to show that uniform contractivity holds. Since in this example the input appears additively, uniform contractivity is simply the requirement that the unforced system (u~0) is contractive. Thus, the main step will be to establish the following technical result, see the Material and Methods: for all t §0, and e T , k 1 , k 2 , and d are arbitrary positive constants, is contracting.
Appealing to Theorem 2, we then have the following immediate Corollary: Proposition 2. For any given nonnegative periodic input u of period T, all solutions of system (12) converge exponentially to a periodic solution of period T.
In the following sections, we introduce a matrix measure that will help establish contractivity, and we prove Proposition 1. We will also discuss several extensions of this result, allowing the consideration of multiple driven subsystems as well as more general nonlinear systems with a similar structure. (A graphical algorithm to prove contraction of generic networks of nonlinear systems can also be found in [18] where this transcriptional module is also studied.)

Proof of Proposition 1
We will use Theorem 2. The Jacobian matrix to be studied is: As matrix measure, we will use the measure m P,1 induced by the vector norm Px j j 1 , where P is a suitable nonsingular matrix. More specifically, we will pick P diagonal: where p 1 and p 2 are two positive numbers to be appropriately chosen depending on the parameters defining the system. Figure 2. A schematic diagram of the transcriptional system modeled in (12). As explained in [29], the transcriptional component takes as input the concentration of protein Z and gives as output the concentration of protein X . The downstream transcriptional module takes as input the concentration of protein X . doi:10.1371/journal.pcbi.1000739.g002 It follows from general facts about matrix measures that where m 1 is the measure associated to the . j j 1 norm and is explicitly given by the following formula: Observe that, if the entries of J are negative, then asking that m 1 (J)v0 amounts to a column diagonal dominance condition.
(The above formula is for real matrices. If complex matrices would be considered, then the term J jj should be replaced by its real part <fJ jj g.) Thus, the first step in computing m P,1 J ð Þ is to calculate PJP {1 : Using (17), we obtain: Note that we are not interested in calculating the exact value for the above measure, but just in ensuring that it is negative. To guarantee that m P,1 J ð Þv0, the following two conditions must hold: Thus, the problem becomes that of checking if there exists an appropriate range of values for p 1 , p 2 that satisfy (20) and (21) simultaneously.
The left hand side of (21) can be written as: which is negative if and only if p 1 vp 2 . In particular, in this case we have: The idea is now to ensure negativity of (20) by using appropriate values for p 1 and p 2 which fulfill the above constraint. Recall that the term e T {y §0 because of the choice of the state space (this quantity represents a concentration). Thus, the left hand side of (20) becomes The next step is to choose appropriately p 2 and p 1 (without violating the constraint p 2 wp 1 ). Imposing p 2 =p 1~1 ze, ew0, (23) becomes Then, we have to choose an appropriate value for e in order to make the above quantity uniformly negative. In particular, (24) is uniformly negative if and only if We can now choose Notice that c 2 depends on both system parameters and on the elements p 1 , p 2 , i.e. it depends on the particular metric chosen to prove contraction. This completes the proof of the Proposition.

Generalizations
In this Section, we discuss various generalizations that use the same proof technique.
Assuming X activation by enzyme kinetics. The previous model assumed that X was created in proportion to the amount of external signal u(t). While this may be a natural assumption if u(t) is a transcription factor that controls the expression of X , a different model applies if, instead, the ''active'' form X is obtained from an ''inactive'' form X 0 , for example through a phosphorylation reaction which is catalyzed by a kinase whose abundance is represented by u(t). Suppose that X can also be constitutively deactivated. Thus, the complete system of reactions consists of X ze 'Y , Figure 3. Entrainment of the transcriptional module (12). Time in minutes on the x-axis. The state of the system (green), y, is entrained to both u(t)~1:5z sin (0:1t) and to a repeating 0,1 f g sequence. System parameters are set to: d~3, k 1 = 1, together with where the forward reaction depends on u. Since the concentrations of X 0 zX zY must remain constant, let us say at a value X tot , we eliminate X 0 and have: We will prove that if u t ð Þ is periodic and positive, i.e. u t ð Þ §u 0 w0, then a globally attracting limit cycle exists. Namely, it will be shown, after having performed a linear coordinate transformation, that there exists a negative matrix measure for the system of interest.
Consider, indeed, the following change of the state variables: x t~x zy: The system dynamics then becomes: As matrix measure, we will now use the measure m ? induced by the vector norm . j j ? . (Notice that this time, the matrix P is the identity matrix).
Given a real matrix J, the matrix measure m ? J ð Þ is explicitly given by the following formula (see e.g. [23]): (Observe that this is a row-dominance condition, in contrast to the dual column-dominance condition used for m 1 .) Differentiation of (28) yields the Jacobian matrix: Thus, it immediately follow from (29) that m ? J ð Þ is negative if and only if: The first inequality is clearly satisfied since by hypotheses both system parameters and the periodic input u t ð Þ are positive. In particular, we have: {u t ð Þ{dz d j jƒ{u 0 :~{c 2 1 ; By using (27) (recall that e T {y §0), the right hand side of the second inequality can be written as: The contraction property for the system is then proved. By means of Theorem 2, we can then conclude that the system can be entrained by any periodic input. Simulation results are presented in Figure 4, where the presence of a stable limit cycle having the same period as u t ð Þ is shown.
Multiple driven systems. We may also treat the case in which the species X regulates multiple downstream transcriptional modules which act independently from each other, as shown in Figure 5. The biochemical parameters defining the different downstream modules may be different from each other, representing a situation in which the transcription factor X regulates different species. After proving a general result on oscillations, and assuming that parameters satisfy the retroactivity estimates discussed in [29], one may in this fashion design a single input-multi output module in which e.g. the outputs are periodic functions with different mean values, settling times, and so forth.
We denote by e 1 , . . . ,e n the various promoters, and use y 1 , . . . ,y n to denote the concentrations of the respective promoters complexed with X . The resulting mathematical model becomes: We consider the corresponding system with no input first, assuming that the states satisfy x(t) §0 and 0ƒy i (t)ƒe T,i for all t,i.
Our generalization can be stated as follows: Proposition 3. System (32) with no input (i.e. u(t)~0) is contracting. Hence, if u(t) is a non-zero periodic input, its solutions exponentially converge towards a periodic orbit of the same period as u(t).
Proof. We only outline the proof, since it is similar to the proof of Proposition 2. We employ the following matrix measure: where P :~p  w0, Vi~1, . . . ,nz1).
In this case, Hence, the nz1 inequalities to be satisfied are: and Clearly, the set of inequalities above admits a solution. Indeed, the left hand side of (38) can be recast as which is negative definite if and only if p 1 =p iz1 v1 for all i~1, . . . ,n. Specifically, in this case we have Clearly, such inequality is satisfied if we choose e 1,iz1 sufficiently small; namely: Þk 2 e T,i : Following a similar derivation to that of the previous Section, we can choose Þk 2 e T,i . In this case, we have: Thus, m J ð Þv{c 2 , where c 2~m in i c i f g, i~1, . . . ,nz1: The second part of the Proposition is then proved by applying Theorem 2. In Figure 6 the behavior of two-driven downstream transcriptional modules is shown. Notice that both the downstream modules are entrained by the periodic input u t ð Þ, but their steady state behavior is different.
Notice that, by the same arguments used above, it can be proven that ð39Þ is contracting. Transcriptional cascades. A cascade of (infinitesimally) contracting systems is also (infinitesimally) contracting [6], [30] (see Materials and Methods for an alternative proof). This implies that any transcriptional cascade, will also give rise to a contracting system, and, in particular, will entrain to periodic inputs. By a transcriptional cascade we mean a system as shown in Figure 7. In this figure, we interpret the intermediate variables  (26). Time in minutes on the x-axis. The system state (green), y, is entrained to the periodic input (blue): u(t)~1:5z sin (0:1t). The zoom on t [ 0,10 ½ min highlights that trajectories starting from different initial conditions converge towards the attracting limit cycle. System parameters are set to: k 1~0 :5, k 2~5 , X tot~1 , e T~1 , d~20. doi:10.1371/journal.pcbi.1000739.g004 ð35Þ ð36Þ Y i as transcription factors, making the simplifying assumption that TF concentration is proportional to active promoter for the corresponding gene. (More complex models, incorporating transcription, translation, and post-translational modifications could themselves, in turn, be modeled as cascades of contracting systems.) More abstract systems. We can extend our results even further, to a larger class of nonlinear systems, as long as the same general structure is present. This can be useful for example to design new synthetic transcription modules or to analyze the entrainment properties of general biological systems. We start with a discussion of a two dimensional system of the form: In molecular biology, a(x) would typically represent a nonlinear degradation, for instance in Michaelis-Menten form, while the function f represents the interaction between x and y. The aim of this Section is to find conditions on the degradation and interaction terms that allow one to show contractivity of the unforced (no input u) system, and hence existence of globally attracting limit cycles. We assume that the state space C is compact (closed and bounded) as well as convex. Since the input appears additively, we must prove contractivity of the unforced system.  Lf Lx .
Notice that the last condition is automatically satisfied if Lf Lx v0, because La Lx w0.
As before, we prove contraction by constructing an appropriate negative measure for the Jacobian of the vector field. In this case, the Jacobian matrix is: Once again, as matrix measure we will use: with and p 1 ,p 2 w0 appropriately chosen. Using (42) we have Following the same steps as the proof of Proposition 1, we have to show that: Clearly, if Lf =Lyw0 for every x,y [ C and p 1 vp 2 , the first inequality is satisfied, with To prove the theorem we need to show that there exists p 1 vp 2 and c 2 2 satisfying (46). For such inequality, since Lf =Lx does not change sign in C by hypothesis, we have two possibilities: In the first case, the right hand side of (46) becomes Choosing p 2 =p 1~1 ze, with ew0, we have: In the second case, the right hand side of (46) becomes Again, by choosing p 2 =p 1~1 ze, with ew0, we have the following Figure 7. Transcriptional cascade discussed in the text. Each box contains the transcriptional module described by (12). doi:10.1371/journal.pcbi.1000739.g007 Figure 6. Entrainment of two-driven transcriptional modules. Time in minutes on the x-axis. Outputs Y 1 (top) and Y 2 (bottom) of two transcriptional modules driven by the external periodic input u(t)~1:5z sin (t). The parameters are set to: d~0:01, k 11~1 0, k 21~1 0, e T,1~1 for module 1 and k 12~0 :1, upper bound for the expression in (48): Thus, it follows that m P,1 J ð Þv{c 2 provided that the above quantity is uniformly negative definite. Since, by hypotheses, The proof of the Theorem is now complete.
From a biological viewpoint, the hardest hypothesis to satisfy in Theorem 3 might be that on the derivatives of f x,y ð Þ. However, it is possible to relax the hypothesis on Lf =Lx if the rate of change of a x ð Þ with respect to x, i.e. La=Lx, is sufficiently larger than Lf =Lx. In particular, the following result can be proved.
Theorem 4. System (40), without inputs u, evolving on a convex compact set, is contractive provided that: Proof. The proof is similar to that of Theorem 3. In particular, we can repeat the same derivation to obtain again inequality (46). Thence, as no hypothesis is made on the sign of Lf =Lx, choosing Thus, it follows that, if La=Lx §2 Lf =Lx j j, then A c 2 such that m P,1 J ð Þv{c 2 , implying contractivity. The above condition is satisfied by hypotheses, hence the theorem is proved.
Remarks. Theorems 3 and 4 show the possibility of designing with high flexibility the self-degradation and interaction functions for an input-output module.
This flexibility can be further increased, for example in the following ways: N Results similar to that of the above Theorems can be derived (and also extended) if some self degradation rate for y is present in (40), i.e.
N Theorem 3 and Theorem 4 can also be extended to the case in which the X -module drives more than one downstream transcriptional modules.

Applications to synthetic biology
We introduced above a methodology for checking if a given transcriptional module can be entrained to some periodic input. The aim of this section is to show that our methodology can serve as an effective tool for designing synthetic biological circuits that are entrained to some desired external input.
In particular, we will consider the synthetic biological oscillator known as the Repressilator [31], for which an additional coupling module has been recently proposed in [32]. A numerical investigation of the synchronization of a network of non-identical Repressilators was independently reported in [33].
We will show that our results can be used to isolate a set of biochemical parameters for which one can guarantee the entrainment to any external periodic signal of this synthetic biological circuit. In what follows, we will use the equations presented in [32] to model the Repressilator and the additional coupling model. Entrainment using an intra-cellular auto-inducer. The Repressilator is a synthetic biological circuit that consists of three genes that inhibit each other in a cyclic way [31]. As shown in Figure 8, gene lacI (associated to the state variable c in the model) expresses protein LacI (C), which inhibits the transcription of gene tetR (a). This translates into protein TetR (A), which inhibits transcription of gene cI (b). Finally, the protein CI (B), translated from cI, inhibits expression of lacI, completing the cycle.
In Figure 9 a modular addition to the three-genes circuit is presented. The module was first presented in [32] and makes the Repressilator circuit sensitive to the concentration of the autoinducer (labeled as S in the model) which is a small molecule that can pass through the cell membrane. Specifically, the module makes use of two proteins: (i) LuxI, which synthesizes the auto-inducer; (ii) LuxR, with which the auto-inducer synthesized by LuxI forms a complex that activates the transcription of various genes.
We model the above circuit with the simplified set of differential equations proposed in [32]. Specifically, the dynamics of the mRNAs are _ a a~{az a Notice that the above equations are dimensionless. This is done by: (i) measuring time in units of mRNA lifetime (which is assumed equal for the three genes), and (ii) expressing the protein levels in units of their Michaelis constant. The parameter a represents the dimensionless transcription rate in the absence of self-repression, while k denotes the maximum contribution of the auto-inducer to the expression of lacI.
The dynamics of the proteins are described by The parameters b A , b B , b C represent the ratios between the mRNAs and the respective proteins' lifetimes and d A , d B , d C represent the protein decay rate. The last differential equation of the model from [32] keeps track of the evolution of the intra-cellular auto-inducer. It is assumed that the proteins TetR and LuxI have equal lifetimes. This in turn implies that the dynamics of such proteins are identical, and hence one uses the same variable to describe both protein concentrations. Thus, the dynamics of the auto-inducer are given by: where k s0 is the rate of degradation of S.
We now model the forcing on the intracellular auto-inducer concentration by adding an external input u t ð Þ to the above dynamical equation. The equation for S becomes: where g can be thought as a diffusion rate.
We will now use the analytical methodology developed in the previous sections, to properly tune the biochemical parameters of the Repressilator circuit, whose mathematical model consists of the set of differential equations (53), (54), (55), so that it shows entrainment to the periodic input u t ð Þ. That is, the measured output (e.g. cI), oscillates asymptotically with a period equal to that of u t ð Þ. Of course, the periodic orbit of the output will depend on the particular choice of the parameters.
In what follows, we assume that all the system parameters can be varied except for the self-degradations that we assume to be fixed as, in practice, they are difficult to modify.
In this case, the Jacobian matrix to be studied is The matrix measure that we will use to prove contraction is where P is a 7|7 diagonal matrix having on the main diagonal the positive arbitrary scalars p i . Computation of PJ JP {1 yields {1z It is easy to check that the nonlinear terms in the above equations satisfy the following inequalities: ð56Þ for all x §0. Hence, the system of inequalities (58a)-(58g) are satisfied, if the following set is fulfilled: {1z {1z The system can then be proved to be contracting for a given set of biochemical parameters, if there exists a set of scalars p i , i~1 . . . 7 satisfying the above inequalities. For example, if the repressilator parameters are chosen so that then it is trivial to prove that, for any constant value p pw0, the set of scalars p i~ p p, for i~1, . . . ,7, satisfies (59a)-(59g). Indeed, in Figure 10 we provide a set of biochemical parameters for which the circuit is contracting and shows entrainment to the periodic input u(t)~1:5z1:5 sin (0:5t). (These parameters, except for the maximal transcription rate a, are in the same ranges as those used in [31], [32]. These parameters are also close to those used in [33] and [34]. The reason for picking an a much smaller than in [32], is that we need to slow down transcription so as to eliminate intrinsic oscillations; otherwise the entrainment effect cannot be shown. This lowering of a by two orders of magnitude is also found in other works, for example in [35], where the same model is studied, with a somewhat larger but of the same order of magnitude as here.) Note that using the set of inequalities (59a)-(59g) as a guideline, it is possible to find other parameter regions where the system is still contracting but exhibit some other desired properties. For instance, to tune (e.g. increase) the amplitude of the output oscillations shown in Figure 10, a possible approach can be that of increasing the biochemical parameter k so as to make stronger the effect of the auto-inducer on the dynamics of the gene cI (variable c(t) in the model).
Again we can prove that the set of inequalities (59a)-(59g) is satisfied for k arbitrarily large, if we set p i~ p p, for i~1, . . . ,6 and choose p 7 such that p p p 7 kv1{M, and {k s0 {gz p 7 p p k s1 ƒ{c 2 7 : Now, due to biochemical constraints the parameter k s1 is considerably smaller than k s0 and g (in our simulations the ratio is of about two orders of magnitude). Therefore, whatever the value of k, it suffices to set p p~1 and p 7~1 0kze, with e being a positive arbitrary constant, to get a solution to (59a)-(59g) and hence guarantee the system to be contracting. Figure 11 shows the behavior of the system output with the modified parameters confirming that with this choice of parameters the oscillation amplitude is indeed larger as expected.
Observe the nonlinear character of the oscillation depicted in Figure 11, which is reflected in the lack of symmetry in the behavior at minima and maxima of cI(t). Our theory predicts the existence (and uniqueness) of such a nonlinear oscillations. None of the usual techniques, based on linear analysis, can explain such behavior.
Entrainment using an extra-cellular auto-inducer. We now consider the case in which the extracellular auto-inducer can change due to an external signal as well as diffusion from intracellular auto-inducer, as represented in Figure 9. A new variable must be introduced, to keep track of the extracellular auto-inducer concentration. The only difference in the new model with respect to the previous one is that the differential equation for S becomes: Notice that the parameter g measures the diffusion rate of the auto-inducer across the cell membrane, i.e. g~sA=V c , with s representing the membrane permeability, A its surface area and V C the cell volume. In the above equation, S e denotes the concentration of the extra-cellular auto-inducer, whose dynamics are given by: where g ext~s A=V ext , with V ext denoting the total extracellular volume, while k se stands for the decay rate.
In analogy with the previous section, we will ensure entrainment of the dynamical system consisting of (53), (54), (61), (62), by tuning the biochemical parameters of this new circuit. Again, the guidelines for engineering the parameters will be provided by the tools developed in the previous sections.
Following the schematic of the previous section, we will prove that there exists c [ R{ 0 f g and a 8|8 constant diagonal matrix P P, such that m P P,? J ð Þƒ{c 2 , where J is the system Jacobian. If we denote with p i , i~1, . . . ,8 the diagonal elements of P P, we obtain the following block-structure for the matrix P PJ P P {1 : Again, we can find sets of biochemical parameters in order to satisfy the above inequalities and hence ensure global entrainment of the circuit to some external input. For example, if we set then, as in the previous section, it is trivial to show that setting all p i to the same identical value satisfies the set of inequality required to prove contraction and hence guarantees entrainment. Notice that the last constraint in (66) is automatically satisfied by the physical (i.e. positivity) constraints on the system parameters.
In Figure 12, the behavior of the circuit is shown with the parameters chosen so as to satisfy the constraints given in (66).
Entraining a population of Repressilators. Consider, now, a population of N Repressilator circuits, which are coupled by means of an auto-inducer molecule. We can think of such a network as having an all-to-all topology, with the coupling given by the concentration of the extracellular auto-inducer, S e . The aim of this section is to show that the methodology proposed here can also be used as an effective tool to guarantee the synchronization of an entire population of biochemical oscillators onto some entraining external periodic input.
We denote with the subscript i the state variables of the i-th circuit in the network, which is modelled using the equations reported in [32] as: : ð67Þ Figure 13 shows a simulation of a population of Repressilators modeled as in (67), with biochemical parameters tuned as in the previous Section: all the circuits composing the network evolve asymptotically towards the same synchronous evolution, which has period equal to that of the input signal u(t). The interested reader is referred to the Materials and Methods for the proof.

Materials and Methods
All simulations are performed in MATLAB (Simulink), Version 7.4, with variable step ODE solver ODE23t. Simulink models are available upon request. The proofs of the results are as follows.

K-reachable sets
We will make use of the following definition: Definition 1. Let Kw0 be any positive real number. A subset C5R n is K-reachable if, for any two points x 0 and y 0 in C there is some continuously differentiable curve c : 0,1 ½ ?C such that: For convex sets C, we may pick c(r)~x 0 zr(y 0 {x 0 ), so c'(r)~y 0 {x 0 and we can take K~1. Thus, convex sets are 1-reachable, and it is easy to show that the converse holds as well.
Notice that a set C is K-reachable for some K if and only if the length of the geodesic (smooth) path (parametrized by arc length), connecting any two points x and y in C, is bounded by some multiple K 0 of the Euclidean norm, y{x j j 2 . Indeed, reparametrizing to a path c defined on 0,1 ½ , we have: Since in finite dimensional spaces all the norms are equivalent, then it is possible to obtain a suitable K for Definition 1. Remark 1. The notion of K-reachable set is weaker than that of convex set. Nonetheless, in Theorem 5, we will prove that trajectories of a smooth system, evolving on a K-reachable set, converge towards each other, even if C is not convex. This additional generality allows one to establish contracting behavior for systems evolving on phase spaces exhibiting ''obstacles'', as are frequently encountered in path-planning problems, for example. A mathematical example of a set with obstacles follows. Example 1. Consider the two dimensional set, C, defined by the following constraints: x 2 zy 2 §1, x §0, y §0 : Clearly, C is a non-convex subset of R 2 . We claim that C is K-reachable, for any positive real number Kw 2 p . Indeed, given any two points a and b in C, there are two possibilities: either the segment connecting a and b is in C, or it intersects the unit circle. In the first case, we can simply pick the segment as a curve (K~1). In the second case, one can consider a straight segment that is modified by taking the shortest perimeter route around the circle; the length of the perimeter path is at most 2 p times the length of the omitted segment. (In order to obtain a differentiable, instead of merely a piecewise-differentiable, path, an arbitrarily small increase in K is needed.)  Proof of Theorem 1 We now prove the main result on contracting systems, i.e. Theorem 1, under the hypotheses that the set C, i.e. the set on which the system evolves, is K-reachable.
Theorem 5. Suppose that C is a K-reachable subset of R n and that f (t,x) is infinitesimally contracting with contraction rate c 2 . Then, for every two solutions x(t)~Q(t,0,j) and z(t)~Q(t,0,f) it holds that: Proof. Given any two points x 0 ð Þ~j and z 0 ð Þ~f in C, pick a smooth curve c : 0,1 ½ ?C, such that c 0 ð Þ~j and c 1 ð Þ~f. Let y t,r ð Þ~Q(t,0,c r ð Þ, that is, the solution of system (1)  where J(y t,r ð Þ,t)~L f Lx (y t,r ð Þ,t). Using Coppel's inequality [28], Hence, we obtain Now, using (70), the above inequality becomes: The Theorem is then proved.
Proof of Theorem 1. The proof follows trivially from Theorem 5, after having noticed that in the convex case, we may assume K~1.

Proof of Theorem 2
In this Section we assume that the vector field f is T-periodic and prove Theorem 2.
Theorem 6. Suppose that: N C is a closed K-reachable subset of R n ; N f is infinitesimally contracting with contraction rate c 2 ; N f is T-periodic; N Ke {c 2 T v1.
Then, there is an unique periodic solution a(t) : ½0,?)?C of (1) having period T. Furthermore, every solution x(t), such that x 0 ð Þ~j [ C, converges to a t ð Þ, i.e. x(t){a(t) j j ?0 as t??. Proof. Observe that P is a contraction with factor Ke {c 2 T v1: P(j){P(f) j j ƒKe {c 2 T j{f j j for all j,f [ C, as a consequence of Theorem 5. The set C is a closed subset of R n and hence complete Figure 11. Increasing the amplitude of oscillations for the model described by (53), (54), (55). Time (minutes) on the x-axis. Behavior of cI when: (i) the input u(t)~0:4z0:4 sin (0:5t) is applied; (ii) no forcing is present. System parameters are the same as that used in Figure 10, except k~15. doi:10.1371/journal.pcbi.1000739.g011 as a metric space with respect to the distance induced by the norm being considered. Thus, by the contraction mapping theorem, there is a (unique) fixed point j j of P. Let a(t) :~Q(t,0, j j). Since a(T)~P( j j)~ j j~a(0), a(t) is a periodic orbit of period T. Moreover, again by Theorem 5, we have that x(t){a(t) j j ƒKe {c 2 t j{ j j ?0. Uniqueness is clear, since two different periodic orbits would be disjoint compact subsets, and hence at positive distance from each other, contradicting convergence. This completes the proof.
Proof of Theorem 2. It will suffice to note that the assumption Ke {c 2 T v1 in Theorem 6 is automatically satisfied when the set C is convex (i.e. K~1) and the system is infinitesimally contracting.
Notice that, even in the non-convex case, the assumption Ke {c 2 T v1 can be ignored, if we are willing to assert only the existence (and global convergence to) a unique periodic orbit, with some period kT for some integer kw1. Indeed, the vector field is also kT-periodic for any integer k. Picking k large enough so that Ke {c 2 kT v1, we have the conclusion that such an orbit exists, applying Theorem 6.

Cascades
In order to show that cascades of contracting systems remain contracting, it is enough to show this, inductively, for a cascade of two systems.
Consider a system of the following form: where x(t) [ C 1 (R n 1 and y(t) [ C 2 (R n 2 for all t (C 1 and C 2 are two K-reachable sets). We write the Jacobian of f with respect to x as A(t,x)~L f Lx (t,x), the Jacobian of g with respect to x as B(t,x,y)~L g Lx (t,x,y), and the Jacobian of g with respect to y as C(t,x,y)~L g Ly (t,x,y), We assume the following: Figure 12. Simulation of the Repressilator forced by some extra-cellular molecule. Time (minutes) on the x-axis. Behavior of cI when the input u(t)~0:4z0:4 sin (0:5t) is applied. Notice that when no forcing is present, the steady state behavior is non-oscillatory. System parameters are: b A~bB~bC~1 , d A~dB~dC~1 :1, a~1:5, k~0:5, k s0~1 , k s1~0 :01. doi:10.1371/journal.pcbi.1000739.g012 Figure 13. Synchronization of Repressilators. Behavior of a population of Repressilator modeled as in (80). Time (minutes) on x-axs. Notice that all the circuits synchronize with a steady-state evolution having the same period as u t ð Þ~0:4z0:4 sin 0:5t ð Þ. System parameters are chosen as in Figure 11, with g ext~0 :1. doi:10.1371/journal.pcbi.1000739.g013 _ y y~g y,v t ð Þ ð Þ, ð74Þ with v t ð Þ being an exogenous input, satisfies the contractivity assumptions of the above Section. Then, consider the interconnection of N identical systems which interact through the variable y as follows: Suppose that ½x 1 (t), . . . ,x N (t),y(t) is a solution of (75) defined for all t §0, for some input u(t). Then, we have the synchronization condition: x i (t){x j (t)?0, as t?z?. Indeed, we only need to observe that every pair ½x i (t),y(t) is a solution of (74) with the same input v(t)~X Furthermore, if u t ð Þ is a T-periodic function, the N interconnected dynamical systems synchronize onto a T-periodic trajectory.
The above principle can be immediately applied to prove that synchronization onto a T-periodic orbit is attained for the Repressilator circuits composing network (67) (see also [19]).
Specifically, let x i :~a i ,b i ,c i ,A i ,B i ,C i ,S i ½ and y~S e ; we have that x 1 , . . . ,x N ,y ½ is a solution of (67). We notice that any pair x i ,y ½ is a solution of the following cascade system ones studied here) should prove invaluable in guiding the search.