# Effect of periodic parametric excitation on an ensemble of force–coupled self–oscillators

## Abstract

We report the synchronization behavior in a one–dimensional chain of identical limit cycle oscillators coupled to a mass–spring load via a force relation. We consider the effect of periodic parametric modulation on the final synchronization states of the system. Two types of external parametric excitations are investigated numerically: periodic modulation of the stiffness of the inertial oscillator and periodic excitation of the frequency of the self–oscillatory element. We show that the synchronization scenarios are ruled not only by the choice of parameters of the excitation force but depend on the initial collective state in the ensemble. We give detailed analysis of entrainment behavior for initially homogeneous and inhomogeneous states. Among other results, we describe a regime of partial synchronization. This regime is characterized by the frequency of collective oscillation being entrained to the stimulation frequency but different from the average individual oscillators frequency.

## 1Introduction

In Nature many life significant processes are regulated by mechanisms of synchronization and entrainment. A properly chosen external stimulus can synchronize an intrinsic biological cycle. A great number of biological processes involve several oscillatory cycles that interact and contribute to an overall system temporal behavior. To give a few examples, the human heart responds sensitively to the forced oscillations modulated by a pacemaker with various frequencies [1], the respiratory rhythm in humans and animals can be entrained to a mechanical ventilator phase [2]. Fibrillar flight muscles of some insects act as a mechanical resonator that is coupled to the insects thorax and wings [3]. A rhythm of flight muscles contractions and wings vibration frequency can be turned by an experimental mechanical or optical stimulation [4]. In all the examples provided so far a self–oscillatory rhythm is modified due to mechanical forces via the coupling to an external mechanical oscillator. A simple theoretical description in terms of interacting nonlinear oscillators would provide a universal approach for exploiting dynamics in this type of systems.

For analysis of time evolution processes simplified models from nonlinear dynamics have been used [6]. In fact, in the framework of nonlinear models one can describe the phenomena that occur in a relatively large class of biological and chemical systems [7]. In this paper we use the analogies found in nature and introduce a non–feedback model of entrainment in a system of coupled nonlinear oscillators. Our model might help to elucidate the basic underlying mechanisms that regulate synchronization of processes in a real biosystem.

Although simplified low–dimensional models can capture inherent mechanisms leading to synchronization, quite often it happens that spatial dynamics plays a crucial role in mediating synchronization [9] in a given system. To unveil the whole variety of possible synchronization scenarios one has to introduce a spatial dimension in the mathematical model of a biosystem.

Collective phenomena in spatially extended systems and their one–dimensional chain analogues have been a subject of a large body of investigations [6]. For instance, models of self–oscillatory systems are widely used to describe biochemical pattern formation processes in spatially extended systems [13]. It is well known that the complex spatiotemporal behavior in a system is largely determined by its dynamical instabilities. There are two basic approaches that are used for numerical and analytical studies of stability in spatially extended systems subjected to external forcing. The first approach implements the mathematical framework provided by the complex Ginzburg–Landau equation (CGLE) [14]. In this case the description of an oscillatory medium is reduced to the study of the amplitude equation with a periodic forcing. Different types of periodic forcing of the CGLE and a rich set of dynamical states and their bifurcations have been studied in [9]. The second approach uses a description of system spatiotemporal dynamics in terms of a large populations of coupled oscillators. Due to a large number of degrees of freedom these systems can display a panoply of dynamical behaviors: from cluster formation [16] to multistability regimes [19] with the coexistence of distinct stable collective modes of oscillations.

Most of the theoretical studies mentioned so far are devoted to investigation of dynamics in globally or locally coupled ensembles of oscillators. The coupling term usually depends on phases or displacements of oscillators. In this work we discuss a new type of coupling via a common force. This force explicitly depends on the displacements as well as on the velocities and the accelerations of individual elements in a large ensemble. The idea to introduce such a coupling has been first proposed in Ref. [20] The model of a large ensemble of self–oscillatory elements coupled via a common force to an inertial oscillator was proposed in [20] as a phenomenological approach to study spontaneous oscillations in an active muscle fiber coupled to mechanical resonator [21]. Each element in the ensemble is described by the Stuart–Landau equation with real coefficient. It is shown that the system undergoes a Hopf bifurcation near to a critical threshold defined in systems parameter space. Various regimes of collective oscillations are reported to exist. In fact, the systems parameter space is divided into several regions with mono– and multi–stable collective modes of oscillatory dynamics. Near to an instability threshold defined in the parameter space there exist completely synchronous, asynchronous and antisynchronous regimes of collective dynamics. In addition, regions with the coexistence of several modes are found.

In this paper we introduce a modification of the model considered in Ref. [20]. We include an external time–dependent stimulation by adding small amplitude time–dependent terms in the parameters expression. We describe the mechanism of entrainment of the frequency of system solutions to the stimulation frequency. The essential feature of the described non–feedback mechanism of entrainment is that the external stimulation enters in the model as a parametric modulation. Two types of parametric excitations are treated: the first type is the periodic modulation of the inertial oscillator stiffness; the second type is the peridic excitation of a natural frequency of a self–oscillatory element. We show that the synchronization scenarios are ruled not only by the choice of parameters of the excitation force but depend on the initial collective state in the ensemble. For both types of the periodic excitations entrainment behavior is studied for homogeneous states. We find sets of frequency–locked solutions that correspond to the resonant driving. We also study the effect of two types of periodic driving on stability of inhomogeneous states. In this work we address several questions: what is the influence of periodic parametric forcing on the stability of a particular collective motion, and how does the features of synchronization compare for different types of parametric excitations in the system.

This paper is organized as follows. In Section 2, we introduce the autonomous model and briefly recall some features of the complex phase dynamics from [20]. Next, in Section 3 we provide the time–dependent model with the first type of parametric modulation and investigate the existence of stable periodic solutions for the case of oscillator in the chain. We explain the numerical procedure that is used for the amplitude–frequency response calculations. Amplitude–frequency responses from numerical simulations and analytical derivation are compared. Next, we present, using numerical simulations, the effect of the periodic forcing on initially inhomogeneous states behavior. We list a variety of regimes of inhomogeneous collective dynamics including partially synchronized states. The second type of parametric modulation is introduced in Section 4. Similarly to the previous section, we provide stability diagrams for the case of oscillatory unit in the chain coupled with an inertial element. We discuss a biological relevance of our findings and give general conclusions in Section 5. In the Appendix an analytical amplitude–frequency relation is derived using a first order quasiperiodic approach .

## 2Model equations

In this section we review the basic model [20] of a chain of active mechanical oscillators coupled to a damped mass–spring oscillator. Each elementary unit in the chain can be modeled by the Stuart–Landau equation. Two variables are defined describing the amplitude and the phase of each th element in a finite chain of elements:

where is the Landau coefficient. The set of equations for describes uncoupled oscillators. In this paper we consider mechanical elements all arranged in a one–dimensional chain and coupled via a common force due to a linear inertial oscillator attached to one end of the chain. In the model [20] a chain of elements has one fixed boundary and is connected to a mass load in a way that the change of the total extension of oscillators is zero. This corresponds to the no flux boundary condition. Consequently, this condition can be recast in the form of the following geometrical constraint imposed on oscillators motion:

Where is the displacement of the th active self–oscillatory element, is a displacement of the linear inertial oscillator. The external force due to the mass–spring load is determined by the equation of motion:

where the parameters and are the mass, characteristic frequency of linear oscillator and the spring damping. The spring stiffness coefficient is . Upon including the mechanical force from (Equation 3) Eqs.(Equation 1) can be written in the complex form [20]:

where and is the scaling parameter. The set of equations (Equation 3)-(Equation 4) represents a system of active oscillators that are coupled via a mean field . Unlike in the Kuramoto model of phase oscillators with global coupling [14] here the collective behavior is governed by a mean field that depends on the velocities and the accelerations of elements in the chain. It was shown in Ref.[20] that oscillatory versus non–oscillatory motions exist in an unfolded space of the system parameters. Namely, oscillatory versus non–oscillatory behavior depends on the values of stiffness, mass of inertial load and the bifurcation parameter . In fact the onset of the oscillatory behavior happens for when the system undergoes a Hopf bifurcation at the critical value of linear frequency . The critical frequency for which transition from a stable non–oscillatory to an oscillating behavior occurs is defined in [20] from the expression: . While oscillatory states can be found for , for higher values of the linear frequency initial oscillations always decay to the stable trivial solution.

A rich bifurcation diagram with various modes of collective oscillatory dynamics exists for near the critical instability threshold [20]. Several regions with homogeneous and inhomogeneous collective modes as well as regions with coexistence of different dynamical regimes are found. Depending on the frequency and one can observe different asymptotic states of system behavior: asynchronous, antisynchronous and synchronous states. For a low value of the stifness the system is set into a synchronized oscillation state. The phase and the amplitude of self–oscillatory elements in a globally synchronized regime are identical. For large positive and for nearly close to the instability threshold value a set of asynchronous states is found. The asynchronous state can be characterized by an inhomogeneous distribution of phases and amplitudes of oscillators along the chain. Remarkably, that while the total chain extension vibrates at one frequency the value of the average oscillators frequency is below it [20]. Localized regions of synchronized motion can be identified along the chain. For large value of the frequency and the system is in an antisynchronous regime. In this state neighboring oscillators move in an antiphase so that the total chain extension remains fixed.

## 3Effect of periodic modulation of stiffness

In this section, we discuss the effect of a non–additive periodic excitation on lock–in solutions stability in the system of self–oscillators coupled to an inertial load. The excitation is taken in the form of periodic modulation of the inertial oscillator stiffness coefficient. First we define a set of equations of motion and then proceed to examining the system stability domains. We study the amplitude dependence for frequencies of excitation near to the one–to–one resonance. Let us introduce a parametric modulation in the stiffness of the inertial oscillator (Equation 3) by using a periodic term with frequency and amplitude . Then the expression for the force is:

where is the time–dependent stiffness coefficient. In our treatment we include the case of a strong magnitude of driving by taking . Upon differentiating both side of Eq. (Equation 4) and substituting the expression for from (Equation 5) the differential equation for the force yields:

Together with Eq. (Equation 4) we now obtain a complete set of differential equations for variables. The displacement of the inertial oscillator is determined from the geometrical constraint (Equation 2). We expect that the stability behavior in system (Equation 4)-(Equation 6) is largely determined by a new type of coupling. As a consequence of the coupling the motion of every self–oscillatory element is affected by a common force. This force is explicitly dependent on time as well as on the oscillators displacement, velocity and acceleration.

### 3.1Entrainment for an initially synchronized regime

Let us first discuss the entrainment behavior for initially homogeneous solutions: all the oscillators in ensemble are moving in phase and with the same frequency of oscillations. The linear stability analysis of the time–independent system in Ref. [20] provides a detailed phase diagram with the regions of stable globally synchronized solutions. Indeed, for the parameters space corresponding to and there exists a large set of stable synchronized solutions of Eqs. (Equation 4)-(Equation 6) for . When the synchronized state is stable it suffices to reduce systems description to the case of oscillator coupled with inertial load. However, when the parametric modulation is introduced in the model the stability problems for the system of and oscillators in the chain are no longer equivalent. The applied periodic excitation can change the stability of the synchronous solution for . Thefore, the system can be driven to a new asymptoticcally inhomogeneous state. In our work we do not analyze the stability of synchronized solutions for the entire parameter space. Instead, we consider entrainment behavior of the synchronized state for a particular choice of parameters: and . We integrate numerically the time-dependent system (Equation 4)-(Equation 6) with and oscillators. It follows that the observed entrainment behavior for the low– and high–dimensional cases matches well for the choice of parameters given above. Therefore, in order to simplify our analytical derivations and to speed up numerical procedure we reduce dynamical description of the system to the low–dimensional case. In the following paragraphs we discuss entrainment behavior for oscillator coupled with a mass load. For initial conditions we chose arbitrary small displacement . Results of numerical simulation are given in Fig. ( ?). We show stability domains in the parameter space defined by the amplitude of driving and the ratio between external frequency and the original limit cycle frequency . Here the stability regions with lock–in solutions are separated from the regions of quasiperiodic and unstable solutions by defined boundaries. Inside each stability region the lock–in solution has rational frequencies ratio . Note that the width of entrainment regions shortens as one proceeds from the resonance to the lower frequency ratios. Several limit cycle solutions are presented in Fig. ( ?) for the resonance and for different amplitudes .

### 3.2Amplitude–frequency response

In this section we present numerical and analytical results on the amplitude–frequency dependence for the system (Equation 4)-(Equation 6). We use numerical analysis of fundamental frequencies [22] to determine the amplitude of quasiperiodic and periodic solutions of the system (Equation 4)-(Equation 6) for the frequencies close to the one–to–one stimulation. For solutions and of (Equation 4)-(Equation 6) we obtain a quasiperiodic approximation near the vicinity of the original limit cycle solution of Eqs. (Equation 4)-(Equation 6) for . By using the substitution of solutions and with their zero time–averaged equivalents: and we write down the quasiperiodic expansions for the new and :

where is a set of time–independent frequencies. We assume that in the above expression the amplitudes and do not dependent on time. In the above expressions the terms are collected according to decreasing order of magnitude. To begin our numerical procedure we introduce a set of new coefficients and that are expressed from the following integrals [22]:

where is the Hanning window. The interval of time is the integration time for every systems trajectory. Next, two sequences of coefficients , and are produced by rearranging and according to the decreasing order of magnitude. In the expansions (Equation 7) we only retain the first two terms that produce reasonably accurate approximation to the solutions and . We calculate the first four largest coefficients of the series and as follows:

Figure ( ?a) illustrates the parametric dependence of the largest coefficient in the expansion ( ?a) on the amplitude and the frequency . Numerical responses are obtained by taking initial conditions for and to be randomly distributed in the vicinity of unperturbed limit cycle for Eqs. (Equation 4)-(Equation 6). For comparison we also calculated the maximum of the numerical solution over the same time interval . Figure ( ?b) shows numerical results for the maximum of the amplitude : . In the plane stability boundaries for the resonance are indicated. The results from two plots are consistent except for the regions near to the stability borders. The amplitude–frequency response is a surface that is smooth inside the borders of stability region and shows discontinuities outside near to the borders of Arnold tongue. Apparently, the discontinuities are due to the emergence of instabilities. Because of the instabilities, the leading frequency components that are defined from the quasiperiodic approximation (Equation 7) in this case might be distinctly different. Consequently, the numerical procedure produces a set of amplitudes which are not smoothly dependent on the external frequency .

We compare numerically estimated amplitude–frequency responses with the results obtained by using quasiperiodic approximation (see Appendix). In Fig. () three numerically calculated amplitude–frequency curves are displayed for different values of the amplitude . To compare these results with the analytical calculations, we display the amplitude–frequency curve obtained from the quasiperiodic approximation of limit cycle solution (see Eq. ( ?) in the Appendix). It is apparent, that for a given choice of parameters the curves from numerical and analytical calculations follow similar nonlinear response behavior inside the stability region (stability boundary is marked by crosses).

### 3.3Entrainment for anti- and asynchronous regimes

In this section we illustrate results of numerical analysis for the periodic forcing of initially spatially inhomogeneous states. The existence and stability of these states for the autonomous system () is controlled by the choice of a mass of inertial oscillator and parameter . For the time–dependent system it is expected that under the influence of a strong driving these states might become unstable and converge to a new asymptotic state. In the limit of a large excitation it can produce significant changes on the regions of stability of existing spatiotemporal regimes. We demonstrate that, meanwhile some collective regimes do not survive new spatially homogeneous states become asymptotically stable. Here we show which conditions on the frequency and the amplitude of external driving have to be satisfied in order for the resulting spatial state to be stable and homogeneous.

In the literature studies of collective dynamics of chains of coupled oscillators are often reduced to their phase dynamics studies [23]. Our system is different from these cases because the amplitude of external force strongly depends on the state of the th oscillator, and therefore, on the collective dynamics in an ensemble. It is not sufficient to confine the investigation only to a phase dynamics study, like in the cases of weakly coupled oscillators [24]. It should be pointed out, that in the limit of a strong excitation the effect on amplitude variation of each individual oscillator has to be taken into account. Namely, for a large relative modulation not only the individual oscillator frequency is affected but also the amplitudes.

First, we proceed by considering the spatiotemporal state of the system (Equation 4)-(Equation 6) with and without external parametric excitation. The time evolution of self–oscillators from an initially inhomogeneous collective state is followed by plotting space–time diagrams. These diagrams display a dynamical collective state of the system. In order to do so, we integrate equations. The initial conditions are prepared by adding random perturbations to the initial zero state. Figure ( ?) shows the space–time diagram of the amplitude of th oscillator on the plane for oscillators. The oscillator number is shown in abscissa and time in ordinate. Figure ( ?a) displays the time evolution of an asynchronous state from the initial time to the final time in the absence of external periodic driving. A closer look reveals that the elements in the half of the ensemble () are grouped in small clusters of oscillators. Each group consists of elements that are phase shifted with respect to the neighboring elements of the same group. As the external stimulation turned on the state in Fig. ( ?a) evolves to a new asynchronous state shown in the panel (b) of the same figure. It presents the time–evolution of an inhomogeneous unstable state for oscillators and the frequency of excitation . It can be easily seen that the distribution of groups of oscillators remains the same as in Fig. ( ?a) but the average self–oscillator frequency is shifted towards the stimulation frequency . The state is not stable, and is sensitive to variations of initial distributions of the positions and the velocities of the oscillators. In fact, the value of numerically estimated Lyapunov exponent is positive in this case. Such regimes are found for the CGLE with an additive periodic forcing [9]. These “turbulent synchronized states” arise in the situation when the forcing is weak and its frequency is chosen nearly the same as a natural oscillation frequency [9]. In this case the dynamical state is characterized by the average frequency of collective excited oscillations being equal to external frequency and its numerical Lyapunov exponent being greater than zero.

Finally, a perfectly synchronized solution is displayed in Fig. ( ?c). Unlike in the previous two plots, here after a short transient time the ensemble is organized into coherent structure with all the individual elements moving in phase and with the frequency . To produce this plot another type of stimulation is used that will be discussed in the following sections.

We will specify the transitions that occur in the system for various stimulation frequencies by plotting the frequency ratio curves versus the stimulation frequency. Figure ( ?) displays two cases for distinct parameter choices. Frequency ratio curves for an initially antisynchronous state and for an initially asynchronous state are shown in Fig. ( ?) (a) and (b). Let us first examine the changes in dynamical behavior if we chose the parameters from the domain that corresponds to an antisynchronous oscillatory regime stable for the time–independent system. This can be easily done by defining the bifurcation parameter and by taking the frequency of the inertial oscillator [20]. We plot in Fig. ( ?a) the frequency ratio curves versus . Several distinct regions that correspond to different collective dynamical states exist. The system remains in the antisynchronous state for an interval . For higher frequencies the solution is driven into the synchronized state. This synchronized state is characterized by all the oscillatory elements moving in–phase and with the same frequency. Upon increasing the excitation frequency the system undergoes a transition to an antisynchronous behavior. We examine now the transitions in frequency that occur for an initially asynchronous state. By adjusting the bifurcation parameter to one insures the existence of an asymptotically stable asynchronous state of the unperturbed system () (Equation 4)-(Equation 6). The frequency of the inertial oscillator is chosen to be the same as defined above. We present results of numerically calculated frequency ratios in Fig. ( ?b). Two curves are presented: the frequency ratio between external driving frequency and the frequency of the total extension versus (squares) and the frequency ratio versus (dotted curve), where is the spatially averaged frequency. One can observe that transitions from one collective state to the other are controlled by changing the frequency . In particular, there exists a plateau that indicates the region of quasiperiodic asynchronous solutions. There are several transient regimes calculated for different values of the external driving : asynchronous states for , the region of coexistence of synchronous and

asynchronous states for at the border with the resonance. A large set of synchronized solutions emerges for the ratio when all the oscillatory elements are phase locked to the frequency ratio regime. Note, that although for the initially antisynchronous state the one–to–one stable entrainment regime is not present (see Fig. ( ?a)) this regime exist for the initially asynchronous state. For a detailed view on this regime we refer to Fig. ( ?) which shows a magnification of the plateau. Numerically calculated nonlinear amplitude response (see Fig. ( ?a)) for a total chain extension and the plateau in Fig. ( ?b) are magnified for ranges of frequencies that lie close to the resonance. In Fig. ( ?c) the evolution of oscillator phases is plotted for . An example of quasiperiodic solution for resonance is presented in Fig. ( ?d). Here, the oscillation of tension force is displayed versus the total chain extension .

Our results suggests, that the final collective behavior depends on the choice of initial collective state of the system. Indeed, as we have seen in Fig. ( ?) the same stimulation frequency results in two distinct final states. One interesting detail, that the stability diagrams indicate the largest set of the entrained solutions. The emergence of these solutions is due to the presence of the multiplicative periodic modulation that produces the harmonic terms in the system equations.

## 4System with parametric modulation of frequency of limit cycle oscillators

In this section we discuss the synchronization dynamics for the system of self–oscillators with a parametric modulation of the frequency of self–oscillatory element. We modify a set of equations (Equation 4)–(Equation 6) by setting and introducing a periodic term with the frequency and the amplitude . The resulting equation for the th oscillator can be written as follows:

Where is the time–dependent frequency modulation term. Consequently, the equation for the force yields:

### 4.1Stability diagram for the synchronized state

Analogously to the previous section where we first discussed the entrainment behavior for the low–dimensional system, here we study first the case of oscillator. By choosing and the linear oscillator frequency below the critical value we assert that the stable globally synchronized state exist for the autonomous system (Equation 8)-(Equation 9) (). We intergate the set of equations (Equation 8)-(Equation 9) with initial conditions picked randomly near to the trivial solution. We carry out numerical stability for a non–zero parametric modulation . In order to reach the limit of a strong excitation we define the amplitude of the modulating term to be comparable with the amplitude of a self–oscillating element frequency . Figure ( ?) shows Arnol’d tongues calculated numerically. Inside the stability regions the frequency ratios are identified for each resonance. Examples of the limit cycle solutions for the unperturbed system () and for the perturbed system with and are shown in Fig. ( ?).

### 4.2Amplitude–frequency response

In Fig. ( ?) the amplitude–response surface is presented from numerical calculations of the maximum of solution for the system (Equation 8)-(Equation 9). The maximum is estimated over the entire interval of integration and obtained from the zero–averaged solution . The borders of stability region are displayed on the plane. The response is smooth inside the one–to–one entrainment stability zone. It strongly shows nonlinear features for large amplitude of the periodic driving: .

### 4.3Entrainment for inhomogeneous states

Now we turn to the case of entrainment behavior for inhomogeneous states. We compute frequency ratio curves for initially inhomogeneous states of a chain of . The stimulation frequency is chosen within the interval in order to include the resonance and the lower ratio resonances. The results are reproduced in Fig. ( ?) for an initially antisynchronous state (a) and for an initially asynchronous state (b). In both cases, one can distinguish two well–pronounced extended plateaus at the and the resonances. The plateaus indicate the perfectly entrained solutions with all the oscillatory units moving in phase and with the frequency of . The remaining solutions that does not belong to the entrainment regions are inhomogeneous states. They display qualitatively different spacial dynamics. Each state is characterized by the individual oscillators forming macroscopic clusters. The motion of the individual units in a cluster is quasiperiodic. One can observe a region of multistability for low values of . Stable coexisting homogeneous and inhomogeneous states can be found inside this region. In the following figure ( ?) the magnification on the one–to–one synchronized state is displayed. This state is obtain from an initially randomly distributed phases of individual elements in the chain. Th amplitude of modulation is and the external frequency . The space–time plot of the synchronized state ca be also viewed in Fig. ( ?c). After a short transient all the oscillators begin to move in phase and with frequency equal to . For the synchronized state shown in Fig. ( ?) the frequency of the total extension and the frequencies of self–oscillatory elements math the stimulation frequency . In Fig. ( ?a) and (b) the amplitude response and the frequency ratio are plotted versus . The periodic orbit in the space of variables and is also shown.

## 5Discussion

In this work we study entrainment behavior in an ensemble of identical limit cycle oscillators coupled via geometrical constraint and via a common force to an inertial linear oscillator. The external force due to a mass–spring load provides a global coupling: it acts on each element of the ensemble. We consider an external parametric modulation as a time–dependent force acting in addition to the applied force due to a mass load The collective dynamics of oscillatory elements undergoes several transitions between different inhomogeneous states. Such transitions are ruled by the frequency and the type of parametric modulations introduced in the model. Based on our description we find families of stable oscillatory synchronized solutions corresponding to the resonant external stimulation of the system. We discuss an influence of two types of parametric excitations: periodic modulation of the stiffness of the inertial load and modulation of the frequency of self–oscillatory element. In particular, we demonstrate that for the spatially synchronized state the second type of excitation leads to larger stability domains than the first type of excitation. Furthermore, we show that with the appropriate parametric modulation introduced an initially inhomogeneous state of the system is driven into the synchronized state with the frequency of collective oscillations equal to the external frequency. Our conclusion is that such globally synchronized behavior can not be achieved from an initially inhomogeneous state by the periodic modulation of the stiffness coefficient. Instead a new type of dynamical state is observed. The emergent spatially inhomogeneous behavior can be described as a partially synchronized state: while every individual element moves with a frequency higher than the frequency of the total chain extension, the sum of displacements of all the elements oscillate at a frequency equal to the external stimulation frequency.

The amplitude–frequency characteristic shows a typical nonlinear response behavior with a fold in the parameter space defined by the driving frequency and the amplitude of the response. Our numerical results show a good agreement with the quasiperiodic approximation inside the one–to–one stable entrainment zone.

The phenomena discussed here have relevance and applications in nature and laboratories. In particular, the autonomous model and its time–dependent analogue provide a phenomenological approach to model oscillatory dynamics of an active muscle fiber. Experimental *in vitro* observations on skinned muscle fiber show evidence of different oscillatory regimes generated by spontaneous contractions of muscle elementary units, sarcomeres [21]. Depending on the conditions of the experiment, oscillatory contractions can undergo a synchronized activity [21] or a spatially inhomogeneous (asynchronous or antisynchronous) activity [26]. In experiments the signalling between adjacent sarcomeres is chemically regulated. The state of activation of one sarcomere is affected by the state of the adjacent sarcomere [21]. Our current model does not consider the effect of coupling between adjacent sites in the chain of elements. In the future we plan to include local coupling between elements and to consider its effect on spatiotemporal dynamics. Detailed study of the corresponding regimes of spatiotemporal behavior is left to future work.

Experimental work and results from theoretical analysis of our phenomenological model promise new interesting insights into collective behavior of oscillations in isolated fibrillar muscle. Periodic modulation of the model parameters considered here might serve as an appropriate setting for modeling active muscle contractions under a specially designed external mechanical and chemical stimulation.

## 6Appendix

In the appendix, we discuss a quasiperiodic approximation [27] used for derivation of the nonlinear amplitude–frequency dependence. The derivation is carried out for the first type of parametric modulation discussed in the paper. It can be seen from Fig. ( ?) that the original unperturbed limit cycle is close to the ellipse centered at zero. Solutions for a non–zero small perturbation are nearly close to the disturbed elliptical form with their centers shifted away from zero. These observations lead us to seek the solutions expansion in form of a finite series of harmonic terms, where the largest contribution terms represent solution for an ellipse. The system (Equation 4)-(Equation 6) can be rewritten in a real form after introducing and variables with a displacement of the oscillator and a time–dependent variable nonlinearly coupled to :

For the non–autonomous system ( ?) a limit cycle solution can be written in the quasiperiodic form [28]:

Where is the frequency of a new limit cycle for the time–dependent system ( ?). The amplitude coefficients in the expansion ( ?) depend on the size of the original unperturbed limit cycle solution and are ordered according to the powers of : are of order , are of order and . If the driving frequency is near to the frequency of the limit cycle solution for the unperturbed system in (Equation 4)-(Equation 6) a reasonable approximation for the lock–in solution can be made by retaining only terms to the accuracy of in Eqs. ( ?). We would not consider the terms of order equal to in the approximation presented here. Coefficients of the expansion ( ?) can be expressed via coefficients in (Equation 7) as follows:

We introduce a small parameter to be a frequency mismatch between the external frequency and the resulting limit cycle frequency. Let us define a linear frequency:

Consequently, one finds that the frequency ratio can be expressed as follows:

In the next step, we use the quasiperiodic approximation ( ?) to derive a final form of the amplitude–frequency response. We assume that the following expansion holds true: , with . For the sake of simplicity, we introduce the rescaling of the parameters:

Upon substituting the quasiperiodic form ( ?) into Eq. (Equation 5) and collecting various terms according to the order of harmonics we obtain a set of nonlinear equations:

where are nonlinear functions of the coefficients in ( ?):

We would not consider the terms of order in the approximation presented here. Let us make the assumption that for the entrained limit cycle solution ( ?) the coefficients in the expansion ( ?) are constant or vary slowly with time. Then one can neglect derivatives from the equations for the coefficients ( ?). By neglecting terms of order one can express the coefficients and :

Let us use the substitution . Our final goal is to obtain a family of amplitude–frequency curves:

After substituting the expressions ( ?) for the coefficients into Eqs. ( ?) and eliminating terms of the order we write down the approximate amplitude–frequency relation:

Finally, what remains to be done is to determine a set of parameters for calculating the relation ( ?). Let us define the parameters in order to compare with the numerical results: . We have used in order to obtain a reasonable fitting to the numerical amplitude–frequency response in Fig. ( ?).

## 7Acknowlegements

The author thanks E. Nicola for many useful comments and suggestions regarding implementation of numerical code. The author is grateful to S. Ares for his comments on the manuscript. The author would like to acknowledge M. Zapotocky for helpful discussions. This research was partially funded by the VolkswagenStiftung Foundation and the Max-Planck-Gesellschaft.

### References

- L. Glass, A. Shrier, in L. Glass, P. Hunter, A. McCulloch,(Eds.),Theory of Heart,Spinger, New York,289 (1991)
- C. Graves, L. Glass, D. Laporta, R. Meloche, A. Grassino, Am. J. Physiol.(Regulat. Integrat. Comp. Physiol. 19) 250 R902 (1986)
- M. H. Dickinson, Current Biology 16 R309 (2006)
- G. Nalbach, Proc. Göttingen Neurobiol. Conf. 16 131 (1988)
- F. -O Lehmann, M. H. Dickinson, J. Exp. Biol. 200 1133 (1997)
- A. Pikovsky, M. Rosenblum, J. Kurths, Synchronization: a Universal Concept in Nonlinear Science,Cambridge University Press, New York,(2001)
- L. Glass, Nature 401 277 (2001)
- H. D. Abarbanel, M. I. Rabinovich, Current Opinion in Neurobiology 11 423 2006
- H. Chate, A. Pikovsky, O. Rudzick, Physica D 131 17 (1999)
- G. V. Osipov, J. Kurths, C. Zhou, Synchronization in Oscillatory Networks, Springer-Verlag, Berlin Heidelberg, (2007)
- V. S. Afraimovich, V. I. Nekorkin, G. V. Osipov, V. D. Shalfeev, Stability, Structures and Chaos in Nonlinear Synchronization Networks, World Scientific, Singapure, (1994)
- S. Rüdiger, E. M. Nicola, J. Casademunt, L. Kramer, Phys. Rep. 447 73 (2007)
- D. Battogtokh, A. Mikhailov, Physica D 90 84 (1996)
- Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer-Verlag, Berlin,(1984)
- P. Coullet, K. Emilsson, Physica D 61 119 (1992)
- S. D. Drendel, N. P. Hors, V. A. Vasiliev, Dynamics of Cell Populations, Nizhny Novgorod University Press, Nizhny Novgorod, (1984)(in Russian)
- T. E. Vadivasova, G. E. Strelkova, V. S. Anishchenko, Phys. Rev. E 63 036225 (2001)
- H. Hakim, W.-J. Rappel, Phys. Rev. A 46 7347 (1992)
- D. Topaj, A. Pikovsky, Physica D 170 118 (2002)
- A. Vilfan, T. Duke, Phys. Rev. Lett. 91 114101 (2003)
- K. Yasuda, Y. Shindo, S. Ishiwata, Biophys. J. 70 1823 (1996)
- , J. Laskar, C. Froeschlé, A. Celletti, Physica D 56 253 (1992)
- H. Daido H, Physica D 91 24 (1996)
- G. B. Ermentrout, N. Kopell, J. Math. Biol. 29 195 (1991)
- N. Fukuda, H. Fujita, T. Fujita and S. Ishiwata, Eur. J. Physiol. 433 1 (1996)
- , N. Okamura, S. Ishiwata, J. Muscle Res.Cell Motil. 9 111 (1988)
- , N. N. Bogoliubov, Y. A. Mitropolsky, Asymptotic Methods in the Theory of Nonlinear Oscillators, Gordon and Breach, New York,(1961)
- K. Tomita, T. Kai, F. Hikami, Prog. Theor. Phys. 57 1159 (1997)