Discrimination of time-dependent inflow properties with a cooperative dynamical system

Many physical, chemical and biological systems exhibit a cooperative or sigmoidal response with respect to the input. In biochemistry, such behavior is called an allosteric effect. Here we demonstrate that a system with such properties can be used to discriminate the amplitude or frequency of an external periodic perturbation or input. Numerical simulations performed for a model sigmoidal kinetics illustrate that there exists a narrow range of frequencies and amplitudes within which the system evolves toward significantly different states. Therefore, observation of system evolution should provide information about the characteristics of the perturbation. The discrimination properties for periodic perturbation are generic. They can be observed in various dynamical systems and for different types of periodic perturbation. \end{abstract}


Introduction
Bistability and hysteresis are commonly observed in physics, chemistry and biology [1][2][3][4][5].Let us assume that a system has two stable states S 1 and S 2 and that an increase in the value of control parameter λ above the threshold λ 1 triggers the transition from S 1 to S 2 , whereas the reverse transition from S 2 to S 1 occurs if the value of the control parameter drops below λ 2 .Such a system can obviously be used as a discriminator of the control parameter value.For example, if the initial state is S 1 and after some time we observe the system in S 2 , then at some point the value of the control parameter necessarily exceeded λ 1 .However, if only time-monotonic changes in the value of the control parameter are considered, then the system discrimination ability is reduced to just two values ; λ 1 and λ 2 .
In this paper we demonstrate the suitability of a dynamical system characterized by sigmoidal kinetics for discrimination-oriented applications, under a new strategy of imposing a periodic perturbation or input on a cooperative system.It has been reported [6][7][8] that periodic perturbations can significantly change the time evolution of a nonlinear system.As a discriminator prototype, we consider a two-variable system in which the inflow of one of the variables is a control parameter.Numerical simulations reveal a non-trivial property of such a system: a marginal change in the inflow parameters (amplitude or frequency) can switch the response of a cooperative system between different branches in the stage diagram.The frequency at which such switching occurs is a monotonic function of the inflow amplitude.Therefore, at a fixed amplitude of periodic inflow, the observation of a transition between different types of oscillatory evolution of the system provides information which allows us to discriminate the inflow frequency.The above discussion does not necessarily limit the range of frequencies that can be discriminated by the observation of transitions between different types of oscillations Similarly, for a fixed frequency of periodic inflow, transitions between different types of system oscillations occur within a narrow range of amplitudes.The transition can be used to discriminate the inflow amplitude, but for the model considered here, the useful range of such discrimination is rather limited.
In numerical simulations, we consider simple system dynamics defined by a single sigmoidal term expressed by a rational function, which is typical for enzymatic reactions [9][10][11][12][13].In such reactions the appearance of sigmoidal kinetic behavior is usually interpreted to be the result of the interaction of substrates with enzymes through positive cooperative binding.Modeling of cooperative binding leads to the Hill equation [9]: where θ is the fraction of ligand binding sites filled, [L] is the ligand concentration, K d is the apparent dissociation constant derived from the mass action law, and n is the Hill coefficient which represents the degree of cooperativity.If n = 1, there is no cooperativity; for n > 1, the cooperativity is positive.Kinetics with sigmoidal behavior are not limited to enzymatic reactions.This also describes the response of various biological systems to external stimuli, including the effect of drug delivery, which is an interesting topic in pharmacology.Among the many experimental studies that have reported sigmoidal behavior, the Hill coefficient n usually has a value between 2 and 4 [9][10][11][12][13][14][15][16][17][18][19].Here we selected n = 3 for the numerical simulations presented below.
The paper is organized as follows.In the next section, we consider a bistable model and study its time evolution as a function of the amplitude and frequency of periodic inflow.We demonstrate how the system can be used as a discriminator and discuss the sources of discrimination errors.In the final section we argue that the observed phenomenon is generic and discuss its potential applications. 2 The response of a model dynamical system to periodic perturbations Let us consider a dynamical system of two variables (x(t), y(t)) defined by a set of differential equations: In Eq.( 2), the last term I(t) = A • (sin(2πf t + φ 0 ) + 1) • Θ(t) describes a periodic inflow of x with frequency f and initial (for t = 0) phase φ 0 .Θ(t) is the Heaviside step function.We assume that the there is no inflow for t < 0, and it is switched on at t = 0.If φ 0 = 3 • π/2, then I(t) is a continuous function.In this case, I(t = 0) = 0.It then increases and finally oscillates.For any other phase, the inflow term is not continuous at t = 0; for example, if φ 0 = π/2 and then I(t = 0) = 2 • A. I(t) then decreases and finally oscillates.The inflow term is always non-negative.For t > 0 the time average of I(t) equals A and is independent of the frequency and the initial phase.If the inflow amplitude A = 0, then (x = 0, y = 0) is the only steady state of Eqs.(2,3) and is stable.In the following analysis, we assume that the stable state of the system without flow is the initial state for the simulated evolution.
Initially, let us consider the time evolution of the system for a constant inflow I(t) = A > 0 for t ≥ 0 (thus, f = 0 and φ 0 = 0).The characteristics of the time evolution depend on the amplitude of the inflow term and on the initial state.In this case, the nullcline g(x, y, t) = 0 is the time-independent line with a definite slope determined by the value of α and a shift which depends on the inflow amplitude A. Figure 1 shows the location of nullclines, calculated for ε = 1, α = 0.55 and a few different values of the inflow.Let us assume that A 1 and A 2 are the amplitudes for which the BN g(x, y) = 0 nullcline is tangential to the sigmoidal-shaped nullcline h(x, y) = 0.The stable stationary states of the system can be located on two branches on the h(x, y) = 0 nullcline.One contains all of the points of the h(x, y) = 0 nullcline located between point (0, 0) and the tangency point (x 1 , y 1 ).We will call it the lower stable branch (LSB).The other is the upper stable branch (USB), and is formed by all of the points of the h(x, y) = 0 nullcline located above (x 2 , y 2 ).The stationary states located on the nullcline between (x 1 , y 1 ) and (x 2 , y 2 ) are unstable.In the case when A < A 2 , the only stationary state is located on the lower stable branch, so the system converges to the stable state y ∞ = lim t→∞ y(t) such that y ∞ ≤ y 1 regardless of the initial state.
Similarly, for A > A 1 , the single BN stationary state is located on the upper stable branch, and for all initial states the system converges to the stable state y ∞ ≥ y 2 .For A 2 ≤ A ≤ A 1 , the stationary state that is approached for t → ∞ depends on the initial state and on the partition of the phase space determined by the separatrices of the saddle point which is located on the middle branch of the h nullcline.This analysis also applies when the frequency of inflow oscillations is very high.In such a case, the flow oscillations are much faster than both the system dynamics and the system responses to the time-averaged value of the inflow A.
For sufficiently slow oscillations of the inflow ( 0 < f ≪ 1), the system can follow the slowly relocating stable state, the position of which varies according to the instantaneous value of the inflow.If the initial state of the system is (x(0) = 0, y(0) = 0) and 2 • A ≤ A 1 , then y(t) ≤ y 1 for all t.Therefore, the system state oscillates along the lower stable branch of the h(x, y) = 0 nullcline with the period defined by the frequency of inflow oscillations.If 2 • A > A 1 , then there are intervals of time within which the system has a single stationary state located on the upper stable branch.
During a single oscillation cycle, there are moments of time t 1 and t 2 at which y(t 1 ) ≤ y 1 and y(t 2 ) ≥ y 2 , and thus oscillations that extend over both stable branches are expected.
For moderate values of f , the system dynamics are too slow to closely follow the changes in the inflow value.In such a case, oscillations around a stable state located in the lower stable branch that extend to the unstable branch, as well as oscillations around a stable state located in the upper stable branch that extend to the unstable branch, should be observed.This is confirmed by numerical simulations.
The complexity of oscillations observed for a constant inflow amplitude A = 0.12 (thus A > A 1 /2 but A < A 1 ) as a function of flow frequency is illustrated in Fig. 2. As discussed above, for the selected amplitude and a low frequency of inflow oscillations, the system dynamics follow the timedependent stationary state and oscillations of y(t) extend over both stable branches of the h(x, y) = 0 nullcline.For intermediate frequencies, oscillations accumulate on the upper stable branch of the nullcline and the minimum value of y(t) increases with frequency.Then, at a certain frequency f c ( for the selected amplitude of oscillations, f c ≈ 0.0312 ), the oscillations switch from the upper to the lower stable branch of the h(x, y) = 0 nullcline.
The transition between oscillations located on different stable branches is quite pronounced and should be easily detected in experiments with a system exhibiting hysteresis.Therefore, it becomes apparent that a cooperative system can discriminate the frequency of a perturbation if its amplitude remains fixed.Numerical simulations have also demonstrated that the frequency of the transition between oscillations on the USB and LSB depends on the phase φ 0 .The right upper corner of Fig. 2 shows two types of oscillations that are observed for f = 0.031.If φ 0 = 0, the system oscillates at the upper stable branch, but if φ 0 = 3π 2 , oscillations around the lower stable branch are seen.Fortunately for the application of this approach to discrimination, the interval of frequencies within which phase-dependent evolution is observed is very narrow.For A = 0.12, it is [0.0306, 0.0312].The width of this interval (∆f ∼ 0.0006) defines the precision in frequency discrimination.
The dynamical system considered here can also be used to discriminate the amplitude of an applied perturbation.Figure 3 shows the time evolution of y(t) for a few values of perturbation amplitude A and a fixed inflow frequency (f = 0.05).As expected, for small amplitudes (A ≤ 0.1122) the oscillations of y(t) are limited to the LSB.For larger amplitudes (0.1122 < A < 0.1361), the range of observed values of y(t) increases, but the oscillations are still anchored on the LSB.Finally, if 0.13617 ≤ A, the oscillations move onto the USB.The transition between the different types of oscillation is quite pronounced and can be used to discriminate the value of amplitude.
Here, similar to the cases illustrated in Fig. 2, we observe a narrow interval of amplitudes ( ∆A ∼ 0.0001) within which the type of oscillation depends on the initial phase.To give a more precise description of system evolution, let us introduce a classification of oscillations based on the minimum and maximum values of y(t) observed over a long time interval for which the evolution has reached a stationary state.We define: and Initially, we used t min = 1000, where t max = t min + 1000.Next we repeated the calculations for t min = 2000.If there is a significant discrepancy in y min , y max obtained for these time intervals, then the procedure is repeated with t min increased by an additional 1000 time units until agreement is attained.
The type of oscillation is classified through the comparison of y min and y max with the values of y 1 and y 2 , as illustrated in Fig. 4. The classification of oscillations is summarized in Table I.
Our discrimination method is based on the determination of conditions in which a small change in f c or A c qualitatively changes the character of the time evolution to force a transition between (I) − 2 and (II) − 1USB type oscillations.Let us assume that we want to measure the unknown inflow frequency and that we can regulate the inflow amplitude.The following procedure can be applied.Initially, we set a low amplitude so the system oscillates on the LSB (type (II) − 1 oscillations).Next, the amplitude is increased up to the moment A z when oscillations of type (II) − 1USB are detected.The frequency of inflow f z can be estimated as f z = G −1 (A z ).This method works for all frequencies greater than f 0 , which corresponds to the tip of the (II) − 1USB region ((f 0 , A 0 )).The accuracy of the estimation depends on the frequency and is high where the amplitude A c is a rapidly increasing function of f c , here for 0.02 ≤ f c ≤ 0.1.This system can also be used to determine the amplitude of inflow when we can control the frequency.Now we set a low frequency and the system exhibits type (III) oscillations.
Next, the frequency is increased up to the moment f y when oscillations of type (I) − 2 are detected.The amplitude of the inflow A y is A y = G(f y ).
Unlike for frequency, the range of discriminated amplitudes does not extend outside the interval [A 0 , A 1 ].
The phase diagrams, similar to that in Fig. 5 but for ε = 1/5 and ε = 5, are shown in Fig. 6 and Fig. 7 respectively.The results are qualitatively identical to those in Fig. 5, suggesting that the described changes in the system oscillations are generic and should also apply to other systems with hysteresis influenced by a periodic perturbation.

Conclusions
We have described how the time evolution of a cooperative system is dependent on the frequency and amplitude of a periodic stimulus.There is a narrow range of these parameters within which the characteristics of this evolution change in a qualitative manner: oscillations around one stable branch change into oscillations on another branch.This phenomenon can be used to determine the amplitude or frequency of an applied perturbation.As for the numerical framework, we evaluated the effect of the rhythmicity of substrate input in a model biochemical system with sigmoidal kinetics, i.This paper describes a system in which the nonlinear term in the kinetic equation for the y(t) variable is described by x 3 1+x 3 term (cf.Eq.( 3)) and the periodic inflow is described by a trigonometric function.We believe that these results are general, and qualitatively similar behavior can be expected in other systems with cooperative characteristics.We performed numerical simulations for a model based on Eqs.( 2,3) but with the inflow term in the form J(t) = A • (tanh (γ • sin(2πf t + φ 0 )) + 1) • Θ(t) for different values of γ.Such periodic inflow becomes a square-like wave for large γ.The phase diagrams that illustrate the type of oscillation as a function of f and A are qualitatively the same, as presented in Figs.(5)(6)(7).We also considered other nonlinear terms in the kinetic equation for y(t), like tanh (x − x 0 ) or 1/(1 + exp (−δ • (x − x 0 )), and obtained similar results.Therefore, we believe that real systems of chemical reactions with hysteresis can be used as discriminators in the manner described above.
The present results can be regarded as a solution to the problem of the optimum stabilization of a system in an unstable state.Let us assume that Eqs.( 2,3) describe the time-dependent progress of a medical treatment where the variable y(t) represents the condition of a patient.The variable x(t) describes the time-dependent concentration of the curing drug.The states on the LSB and USB correspond to an ill and healthy patient, respectively.This simple model seems to realistically describe the basic features of drug therapy.
It predicts that if the inflow of the drug is small, then the patient remains ill.Only a dose higher than a critical dose allows for successful treatment.
However, some drugs are toxic ( such as those used in chemotherapy) and the total dose should be as small as possible.An analysis of the dynamical system presented in Fig. 5 can provide a solution: if we consider the periodic inflow of a drug in the form of Eq.( 1), then the minimum amount of drug required to stabilize the patient in a healthy state corresponds to the bottom corner of the type (II) − 1USB oscillation region -here A 0 ∼ = 0.1 and f 0 ∼ = 0.015.

Figure 3 :
Figure 3: Time evolution of the dynamical system described by Eqs.(2,3) as a function of the inflow amplitude A for a fixed frequency f = 0.05.Tics and numbers on the frequency scale mark transitions between different types of oscillation.The initial phase is φ 0 = 3π/2 for all cases except at the top in the right column, for which φ 0 = 0.The horizontal dashed lines mark the values of y 1 and y 2 .

1LSBFigure 4 :
Figure 4: Geometrical illustration of types of oscillation in the classification based on y min and y max .The dotted lines mark y 1 = 0.09607 and y 2 = 0.65338.

Figure 5
Figure 5 illustrates the regions of parameters (f, A) for which a given oscillation pattern is observed in a model characterized by α = 0.55 and ε = 1.The thick line separates the region of the phase space (f, A) in which class (I)-2 oscillations are observed, from the region where class (II)-1 USB oscillations appear.Let us denote points on this line as (f c , A c ).The amplitude A c , when treated as a function of f c , is a continuous, monotonically increasing function A c = G(f c ).Therefore, the inverse function f c = G −1 (A c ) exists.

1 A 1 / 2 Figure 5 :A 1 / 2 . 1 A 1 / 2 Figure 6 :
Figure 5: Phase diagram showing the type of oscillation as a function of inflow parameters (f, A).The horizontal lines indicate A 1 = 0.16445 and A 1 /2.The model parameters are ε = 1 and α = 0.55 .The thick solid line marks the boundary between oscillation classes (I) − 2 and (II) − 1USB.The transition between these oscillations is used to determine the parameters of inflow.

Figure 7 :
Figure 7: Phase diagram showing the type of oscillation as a function of inflow parameters (f, A).The horizontal lines indicate A 1 = 0.16445 and A 1 /2 .The model parameters are ε = 5 and α = 0.55 .
e. n = 3 in the Hill equation.Using numerical simulations, we separated the phase space of inflow parameters (amplitude and frequency) into regions where specific types of oscillation are observed.The boundary line separating oscillations with significantly different behaviors (type (I) − 2 and type (II) − 1USB oscillations) was identified.The frequency that causes a transition appears in a monotonic function of the inflow amplitude.The system can be used to determine the inflow frequency if we can control the inflow amplitude.It can also be used to determine the inflow amplitude when we can control the frequency.In other words, sigmoidal kinetics with the Hill equation can act as an inflow discriminator.