# Optimal Subharmonic Entrainment

###### Abstract

For many natural and engineered systems, a central function or design goal is the synchronization of one or more rhythmic or oscillating processes to an external forcing signal, which may be periodic on a different time-scale from the actuated process. Such subharmonic synchrony, which is dynamically established when control cycles occur for every cycles of a forced oscillator, is referred to as : entrainment. In many applications, entrainment must be established in an optimal manner, for example by minimizing control energy or the transient time to phase locking. We present a theory for deriving inputs that establish subharmonic : entrainment of general nonlinear oscillators, or of collections of rhythmic dynamical units, while optimizing such objectives. Ordinary differential equation models of oscillating systems are reduced to phase variable representations, each of which consists of a natural frequency and phase response curve. Formal averaging and the calculus of variations are then applied to such reduced models in order to derive optimal subharmonic entrainment waveforms. The optimal entrainment of a canonical model for a spiking neuron is used to illustrate this approach, which is readily extended to arbitrary oscillating systems.

Key words. entrainment, synchronization, oscillators, optimal, time-scales

## 1 Introduction

The synchronization of interacting cyclical processes that evolve on different time-scales plays a fundamental role in many natural phenomena and engineered structures [blekhman88, pikovsky01]. The concept of synchronization is particularly significant to the study of biological systems [strogatz01, granada09enz], which exhibit endogenous oscillations with periods ranging from milliseconds, such as in spiking neurons [izhikevich07], to years, as in hibernation cycles [mrosovsky80]. All life on earth consists of rhythmic systems that are affected by other cyclical processes, such as the daily light and temperature changes that actuate circadian pacemakers [doyle06], and the complex environmental cycles that drive plant growth [mcclung11]. Biological systems may also interact in networks of responsive dynamical units, such as the neurons that constitute circadian oscillators [gonze05] and central pattern generators [coombes05] in the brain, or communicating insects [hartbauer05], for which even the simplest models can produce very complex dynamics [strogatz00].

The process of entrainment, which refers to the dynamic synchronization of an oscillating system to a periodic input, is significant in biology [hanson78, ermentrout84, glass88, granada09enz], with particular relevance in neuroscience [demir97, berke04, sirota08], and is also observed in reactive chemical systems [kuramoto84, aronson86, lev89]. The notion of entrainment is paramount for understanding rhythmic systems, as well as for controlling such systems in an optimal manner [granada09fast, harada10]. Optimal entrainment also has compelling applications in clinical medicine, such as protocols for coping with jet lag [vosko10, sthilaire12], clinical treatments for neurological disorders including epilepsy [kiss08, good09], Parkinson’s disease [hofmann11], and tinnitus [strauss05], and optimization of cardiac pacemakers [naqvi12]. Techniques for controlling the entrainment process can also be used in the design of vibrating mechanical structures [blekhman88, zalalutdinov03] and nanoscale electromechanical devices [feng08, barnes11] that require frequency control or phase locking, and can enable transformational technologies such as neurocomputers [hoppensteadt00] and chaos communication [fischer00]. Various phenomena such as noise-induced synchronization [wacker11], time-scales in synchronization and network dynamics [chen09, chavez05], and transient phenomena [granada09fast, zlotnik13prl] have been examined.

Nonlinear oscillating systems are often studied by transforming the complex dynamic equations that describe their behavior into phase coordinates [ermentrout96, strogatz00], which can also be experimentally established for a physical system when the dynamics are unknown [galan05, tokuda07]. Such models have been studied extensively, with a particular focus on neural [ermentrout96, hoppensteadt99] and electrochemical [kuramoto84, kiss02, nakata09] systems. The reduction of a system model from a complicated set of differential equations to a simple scalar phase coordinate representation is especially compelling from a control-theoretic perspective because it enables a corresponding reduction in the complexity of optimal control problems involving that system. Optimal control of phase models has been investigated with various objectives, such as to alter the spiking of a single neuron using minimum energy inputs [moehlis06] with constrained amplitude [dasanayake11, dasanayake11cdc] and charge balancing [danzl10, dasanayake11cb], as well as to control a network of globally coupled neurons [nabi11]. Several studies have focused on optimal waveforms for entrainment using basic models [schaus06, harada10], and our recent work has resulted in optimal entrainment controls for general nonlinear oscillators that require no knowledge about the initial state or phase of the system [zlotnik11], and can account for uncertainty in oscillation frequency [zlotnik12jne]. These investigations have demonstrated that phase coordinate reduction provides a practical approach to the optimal control of complex oscillating systems.

Previous work on optimal control of the entrainment process has focused on the harmonic case, which corresponds to a one-to-one (1:1) relationship between the frequencies of the stimulus and oscillator. Many physical processes, however, undergo subharmonic : entrainment, which transpires when cycles of the stimulus occur for every cycles of the oscillator [glass88]. Originally examined in the context of loudspeaker dynamics [cunningham51], subharmonic synchronization can emerge among weakly coupled oscillators [ermentrout81, guevara82], and can be induced in forced or injection-locked oscillators to produce entrainment [daryoush89, storti88]. Subharmonic locking phenomena are of interest in a wide range of fields, and neuroscience in particular. Applications exist in magnetoencephalography [tass98], the study of brain connectivity [honey07], dynamic neural regulation [hunter03], as well as in the clinical treatment of epilepsy [hofmann11, kiss08, kriellaars94]. Subharmonic entrainment plays a central role in our understanding of human perception of beat and meter [large96, clayton05, nozaradan11], as well as sound in general, and an ability to affect this phenomenon will lead to innovative therapies for tinnitus [strauss05, tass12]. Indeed, the functional connectivity of the cerebral cortex may be shaped by mutual entrainment of bursting neurons across multiple time scales in a coevolutionary manner [honey07, lajoie11, gutierrez11]. Previous studies have found that subharmonic synchronization phenomena are ubiquitous in biological systems. In fact, respiration and heartbeat in human beings is typically entrained at a 1:4 ratio [schafer98], and evidence exists that human sleep latency is entrained by the lunar cycle [cajochen13, foster08], which is a 1:28 ratio. Other investigations have focused on engineering subharmonic locking in electronic circuits [maffezzoni10, takano07], antenna systems [zarroug95], and voice coil audio systems [bolanos05, noris12].

In this paper, we develop a method for engineering weak, periodic signals that achieve subharmonic entrainment in nonlinear oscillating systems without the use of state feedback. We apply the methods of phase model reduction, formal averaging, and the calculus of variations, which have been used in our previous studies on harmonic entrainment [zlotnik11, zlotnik12jne, zlotnik13prl]. In addition to yielding optimal waveforms for entrainment using weak forcing, this approach allows us to approximate the entrainment regions called Arnold tongues whereby the frequency-locking characteristics of the controlled system are visualized [ermentrout81, schaus06]. We have previously used such graphs to characterize the performance of optimal controls derived using the phase response curve (PRC) of an oscillator when it is applied for harmonic entrainment of the original oscillator in state space [zlotnik12jne]. This crucial validation is extended to the methodology that we apply here to the subharmonic case.

In the following section, we discuss the phase coordinate transformation for a nonlinear oscillator and various methods for computing the PRC. In Section LABEL:secent, we describe how averaging theory is used to study the asymptotic behavior of an oscillating system under subharmonic rhythmic forcing, and in Section LABEL:secminpow we use the calculus of variations to derive the minimum-energy subharmonic entrainment control for a single oscillator with arbitrary PRC. In Section LABEL:secfast, a similar approach is applied to derive the control of fixed energy that produces the fastest subharmonic entrainment of a single oscillator. In Section LABEL:secnmens, we formulate a result on minimum-energy subharmonic entrainment of ensembles of structurally similar nonlinear oscillators, and in Section LABEL:secmaxr we study the dual objective of entraining the largest collection of such oscillators with a control of fixed energy. This is followed by Section LABEL:sechh, where we examine the performance of minimum energy subharmonic entrainment waveforms computed for the Hodgkin-Huxley model [hodgkin52], as well as Section LABEL:secfastsim, where the convergence rate of the system subject to inputs for fast subharmonic entrainment is examined. Finally, in Section LABEL:secdisc we discuss several details and implications of this paper. Throughout the manuscript, important concepts are described graphically using illustrations as well as examples using the phase model of the Hodgkin-Huxley system, which is described in Appendix LABEL:aphh, as a canonical nonlinear oscillator.

## 2 Phase models

The phase coordinate transformation is a model reduction technique that is widely used for studying oscillating systems characterized by complex nonlinear dynamics [kuramoto84], and can also be used for system identification when the dynamics are unknown [kiss05]. Consider a full state-space model of an oscillating system, described by a smooth ordinary differential equation system

\hb@xt@.01(2.1) |

where is the state and is a control. Furthermore, we require that (LABEL:sys1) has an attractive, non-constant limit cycle , satisfying , on the periodic orbit . In order to study the behavior of this system, we reduce it to a scalar equation

\hb@xt@.01(2.2) |

which is called a phase model, where is the phase response curve (PRC) and is the phase associated to the isochron on which is located. The isochron is the manifold in on which all points have asymptotic phase [brown04]. It is standard practice to define (mod ) when the first variable in the state vector attains its maximum over the orbit . This is due to the significant role of mathematical neuroscience in the development of phase model theory. For models of neural oscillators, the first state variable often denotes the membrane potential, which exhibits spiking or relaxation behavior, so that for occur concurrently with successive spikes. The conditions for validity and accuracy of phase reduced models have been determined [efimov10, efimov11], and the reduction is accomplished through the well-studied process of phase coordinate transformation [efimov09], which is based on Floquet theory [perko90, kelley04]. The model is assumed valid for inputs such that the solution to (LABEL:sys1) remains within a neighborhood of .

To compute the PRC, the period and the limit cycle must be approximated to a high degree of accuracy. This can be done using a method for determining the steady-state response of nonlinear oscillators [aprille72] based on perturbation theory [khalil02] and gradient optimization [peressini00]. The PRC can then be computed by integrating the adjoint of the linearization of (LABEL:sys1) [ermentrout96], or by using a more efficient and numerically stable spectral method developed more recently [govaerts06]. A software package called XPPAUT [ermentrout02] is commonly used by researchers to compute the PRC. We employ a technique derived from the method of Malkin [malkin49] in order to compute PRCs, for which details are given in Appendix LABEL:cprc. The PRC of the Hodgkin-Huxley system with nominal parameters, along with the limit cycle, is shown in Figure LABEL:s1f1_prc.

## 3 Fundamental theory of subharmonic entrainment by weak forcing

An essential objective in all entrainment applications is to force the frequency of an oscillator to a desired value. While this can be accomplished using any sufficiently powerful rhythmic signal, it is often desirable to do so using an input that consumes minimum energy, or satisfies another optimization objective. The harmonic (:) case constitutes the canonical entrainment problem, which was examined for arbitrary nonlinear oscillating systems in our previous work [zlotnik11, zlotnik12jne]. The theory of subharmonic (:) entrainment that is presented here is a nontrivial extension of those results.

Our goal is to entrain the system (LABEL:sys2) to a target frequency using a periodic forcing control of frequency , such that cycles of the oscillator occur for every cycles of the input. When such : entrainment occurs, then the target and forcing frequencies satisfy , so that the control input has the form , where is -periodic. From here on, it is assumed that and are coprime integers. In addition, we adopt the weak forcing assumption, i.e., where has unit energy and , so that given this control the state of the original system (LABEL:sys1) is guaranteed to remain in a neighborhood of in which the phase model (LABEL:sys2) remains valid [efimov10]. We then define a slow phase variable by , and call the difference between the natural and target frequencies the frequency detuning. The dynamic equation for the slow phase is then

\hb@xt@.01(3.1) |

where is called the phase drift. In order to study the asymptotic behavior of (LABEL:sys3) it is necessary to eliminate the explicit dependence on time on the right hand side, which can be accomplished by using formal averaging [kuramoto84]. If is the set of -periodic functions on , we can define an averaging operator by

\hb@xt@.01(3.2) |

In addition, let us define the forcing phase and a change of variables . Then the weak ergodic theorem for measure-preserving dynamical systems on the torus [kornfeld82] implies that for any periodic function , the interaction function

\hb@xt@.01(3.3) |

exists as a continuous, -periodic function in . In addition, because both and are -periodic, can be expressed by integrating with respect to or to to yield two equivalent expressions given by

\hb@xt@.01(3.4) | ||||

\hb@xt@.01(3.5) |

In particular, the expression (LABEL:lnm2a) can be written as

\hb@xt@.01(3.6) |

where we define the function

\hb@xt@.01(3.7) |

We henceforth write . At this point, let us establish several important expressions that will be used throughout the following sections. First, we define a function as the : interaction function of with itself by

\hb@xt@.01(3.8) |

By defining an inner product by , the Cauchy-Schwarz inequality yields , and the periodicity of results in . We can then define

\hb@xt@.01(3.9) |

which inherits the properties for all and from the properties of . We will subsequently write and . The expression (LABEL:vnmdef) is important because using in (LABEL:lnm3a) yields

\hb@xt@.01(3.10) |

In addition, using (LABEL:lnm2b) we see that the energy of the function is given by

\hb@xt@.01(3.11) |

The functions , , and , as defined in (LABEL:ynmdf), (LABEL:qdef), and (LABEL:vnmdef), respectively, will appear repeatedly in the subsequent derivations of optimal subharmonic entrainment controls.

As in the case of 1:1 entrainment, the formal averaging theorem [hoppensteadt97] permits us to approximate (LABEL:sys3) by the averaged system

\hb@xt@.01(3.12) |

in the sense that there exists a change of variables that maps solutions of (LABEL:sys3) to those of (LABEL:sys4). A detailed derivation for the : case is given in Appendix B of [zlotnik12jne], and this can be easily extended to the : case. Therefore the weak forcing assumption with allows us to approximate the phase drift equation by

\hb@xt@.01(3.13) |

The averaged equation (LABEL:sys5) is autonomous, and approximately characterizes the asymptotic behavior of the system (LABEL:sys2) under periodic forcing. Specifically, we say that the system is entrained by a control when the phase drift equation (LABEL:sys5) satisfies , which will occur as if there exists a phase that satisfies . When both the control waveform and PRC are non-zero, the function is not identically zero, so when the system is entrained there exists at least one that is an attractive fixed point of (LABEL:sys5). The stable fixed points of (LABEL:sys5), which are the roots of the equation , determine the average phase shift, relative to , at which the oscillation stabilizes from a given initial phase. In addition, we define the phases and at which the interaction function achieves its maximum and minimum values, respectively. In order for entrainment to occur, must hold, so that at least one stable fixed point of exists. Thus the range of the interaction function determines which values of the frequency detuning yield phase locking. These properties are illustrated in Figure LABEL:s3f1_inter1.

Moreover, the interaction function can be used to estimate the values of the minimum root mean square (RMS) energy that results in locking of an oscillator to a given frequency at a subharmonic : ratio using the waveform . This is accomplished by substituting the expression and the relation into the equation and simplifying to obtain

\hb@xt@.01(3.14) |

where is a unit energy normalization of . This equation is then solved for at and to produce linear estimates of boundaries for the regions of pairs that yield entrainment. These regions are known as Arnold tongues, so named after mathematician who first described a similar phenomenon for recurrent maps on the circle (Section 12 of [arnold61]). The RMS energy is used because the boundary of the entrainment region is approximately linear for weak forcing, and yields a clear visualization [ermentrout81, schaus06]. The Arnold tongue boundary estimates obtained using (LABEL:arnold_tg1) can be classified into three different cases that depend on the signs of and , which are listed in Table LABEL:s3t1_arnold and illustrated in Figure LABEL:s3f2_arnold.

Based on the theoretical foundation and fundamental notations presented in this section, we proceed to formulate and solve several design and optimization problems for subharmonic entrainment of one or more oscillating systems. In the following section, we address the canonical problem of establishing subharmonic resonance of a single oscillator to a periodic input of minimum energy at a desired frequency.

## 4 Minimum energy subharmonic entrainment of an oscillator

In many applications described in Section LABEL:secintro, it is desirable to achieve entrainment of an oscillator by using a control of minimum energy. This problem can be formulated as a variational optimization problem where the objective function to be minimized is the control energy , and the design constraint is if and if . This inequality is active when optimal entrainment occurs, and hence can be expressed as the equality constraint

\hb@xt@.01(4.1) | |||||

\hb@xt@.01(4.2) |

We formulate the problem for to obtain the minimum energy control for frequency increase using the calculus of variations [gelfand00]. The derivation of the case where is similar, and results in the symmetric control . The constraint (LABEL:const1a) can be adjoined to the cost using a multiplier , leading to the objective

\hb@xt@.01(4.3) | ||||

where the expression (LABEL:lnm2b) is substituted for . Applying the Euler-Lagrange equation [gelfand00], we obtain the necessary condition for a candidate optimal solution

\hb@xt@.01(4.4) |

Recalling (LABEL:lnmv1), we obtain

\hb@xt@.01(4.5) |

Therefore the constraint (LABEL:const1a) results in

\hb@xt@.01(4.6) |

so the multiplier is given by . Consequently, we can express the minimum-energy control that entrains the system (LABEL:sys2) to a target frequency using a periodic forcing control , where is -periodic and , as

\hb@xt@.01(4.7) |

where is the forcing phase. In practice, we may omit the phase ambiguity or in (LABEL:solnm1) because entrainment is asymptotic. The optimal waveforms for minimum-energy entrainment of the Hodgkin-Huxley system are shown in Figure LABEL:s4f1_mpk for values of , and the corresponding interaction functions are shown in Figure LABEL:s4f1_mpif. Observe that in Figure LABEL:s4f1_mpk, the : minimum-energy control will repeat the : optimal waveform times during the control cycle. As the ratio grows large, the controls converge to a constant given by

\hb@xt@.01(4.8) |

as seen by substituting the limiting expressions as for and from equations (LABEL:ynmdf) and (LABEL:vnmdef) into the solution (LABEL:solnm1).

Finally, by invoking (LABEL:lnmv2), we see that the minimum energy waveform (LABEL:solnm1) for subharmonic entrainment of a single oscillator has energy given by

\hb@xt@.01(4.9) |

We have shown that the minimum energy periodic control that achieves subharmonic entrainment of a single oscillator with natural frequency to a target frequency is a re-scaling of the function given in (LABEL:ynmdf), where is the forcing phase. In the case that , these results reduce to the solution in the harmonic (1:1) case, which is a re-scaling of the PRC [zlotnik11, zlotnik12jne]. In addition, it is important to note that very similar optimal inputs for altering the frequency of oscillating neural systems in the phase model representation have been obtained using other methods in the harmonic (1:1) case [moehlis06, dasanayake11]. The fundamental conclusion is that when the desired change in the frequency of the oscillator is small and the applied input to the system is weak, the optimal control in the harmonic case is a re-scaling of the PRC. Therefore we expect that applying alternative methods [moehlis06, dasanayake11] to compute inputs for minimum energy subharmonic control of oscillators will also result in solutions similar to (LABEL:solnm1).

## 5 Fast subharmonic entrainment of an oscillator

An alternative objective to minimizing control energy is to entrain a system to a desired frequency as quickly as possible using a control of given energy. This problem has been examined for the harmonic (:) case for arbitrary nonlinear oscillating systems in our previous work [zlotnik13prl], and the solution for the subharmonic (:) case extends those results by applying the techniques derived in Section LABEL:secent.

Our goal here is to entrain the system (LABEL:sys2) to a target frequency as quickly as possible by using a periodic control of fixed energy and forcing frequency that satisfies . Employing averaging theory as in Section LABEL:secent yields the phase drift equation (LABEL:sys5), where the interaction function would ideally be of a piecewise-constant form, so that the averaged slow phase converges to a fixed point at a uniform rate from any initial value. However, a discontinuity at would result in an unbounded control , as explained in Lemma 2 of Appendix LABEL:apnum, which makes such a control infeasible in practice. An alternative is to maximize , the rate of convergence of the averaged slow phase in the neighborhood of its attractive fixed point . The calculus of variations can then be used to obtain a smooth optimal candidate solution that also performs well in practice. When the system (LABEL:sys5) is entrained by a control , there exists an attractive fixed point satisfying and , as seen in Figure LABEL:s5f0_inter1. Observe that by inspecting (LABEL:lnm2b), one can write

\hb@xt@.01(5.1) |

where is as defined in (LABEL:ynmdf) and is its derivative, given by

\hb@xt@.01(5.2) |

with . As was done for in Section LABEL:secent, we can define the interaction function of with itself by

\hb@xt@.01(5.3) |

which is maximized at with the maximum value . The periodicity of implies that for all , and . We then define

\hb@xt@.01(5.4) |

which inherits the properties for all and from the function . We will use the notation and . By substituting for and for in (LABEL:lnmv1), we can obtain

\hb@xt@.01(5.5) |

In addition, by combining (LABEL:dlnmdf) and (LABEL:dlnmv1), the energy of the function is given by

\hb@xt@.01(5.6) |

The functions , and , as defined in (LABEL:rrdef) and (LABEL:wwnmdef), respectively, will appear repeatedly in the following derivation of fast subharmonic entrainment controls.

In order to maximize the rate of entrainment in a neighborhood of using a control of energy , the value of should be maximized for values of near , which occurs when is large, as illustrated in Figure LABEL:s5f0_inter1. This results in the following optimal control problem formulation for fast subharmonic entrainment:

\hb@xt@.01(5.7) | ||||

\hb@xt@.01(5.8) | ||||

\hb@xt@.01(5.9) |

The constraints can be adjoined to the objective using multipliers and to yield

\hb@xt@.01(5.10) |

The associated Euler-Lagrange equation is

\hb@xt@.01(5.11) |

and solving for yields the candidate solution

\hb@xt@.01(5.12) |

The multipliers and can be found by substituting (LABEL:sol:cand) into the constraints (LABEL:eq:con1) and (LABEL:eq:con2). This yields the equations

\hb@xt@.01(5.13) | ||||

\hb@xt@.01(5.14) |

where averaging is done with respect to the variable . Because is -periodic, then is as well, as are and in both arguments. Thus one can show, e.g., using Fourier series, that , and also, so that (LABEL:const:2) easily yields

\hb@xt@.01(5.15) |

where can be seen from (LABEL:lnmv2). Substituting this result into (LABEL:const:1) leads to a quadratic equation for given by

\hb@xt@.01(5.16) |

Substituting (LABEL:sol:cand) into the equation (LABEL:dlnmdf), and recalling that , yields

\hb@xt@.01(5.17) |

In particular, , so we choose when solving (LABEL:lfasteq) for in order for the expression (LABEL:fenm4) to be negative in order for the objective in (LABEL:prob:1) to be maximized. It follows that the optimal waveform and multiplier can be obtained from (LABEL:sol:cand), (LABEL:mufasteq) and (LABEL:lfasteq) as

\hb@xt@.01(5.18) |

where phase-locking occurs fastest when the oscillator is in the neighborhood of the phase at the start of entrainment with . For zero frequency detuning, the optimal waveform is a re-scaling of , which is a sum of shifted derivatives of the PRC function. As increases, continuously transforms towards a rescaling of , which is the minimum energy waveform for subharmonic entrainment, as derived in the previous section. This transition reflects the conceptual trade-off between the fast entrainment objective (LABEL:prob:1) and frequency control constraint (LABEL:eq:con2), which can be satisfied only when , as shown in (LABEL:minpvpow). When , these results reduce to the harmonic (1:1) case [zlotnik13prl]. Figure LABEL:s5f1_fk shows the fast subharmonic entrainment controls for the Hodgkin-Huxley model for values of at detunings between and , and Figure LABEL:s5f2_fif shows several corresponding interaction functions. It is important to note that the choice of control energy significantly impacts the control waveform in the case of non-zero detuning, as seen in the panels on the diagonal of Figure LABEL:s5f1_fk.

## 6 Minimum energy subharmonic entrainment of oscillator ensembles

In practice, biological systems exhibit variation in parameters that characterize the system dynamics, which must be taken into account when designing optimal entrainment controls. We approach this issue by modeling an ensemble of systems as a collection of phase models with a common PRC that is derived using a nominal parameter set, and the frequencies span the range of natural frequencies resulting from phase model reduction of systems with parameters in a specified range. A justification of this approach and a sensitivity analysis is provided in Section 5 in [zlotnik12jne]. The following extension of this modeling and control technique to subharmonic (:) entrainment parallels our previous work [zlotnik12jne] by incorporating the theory derived in Section LABEL:secent, and contains a more rigorous optimality proof and additional generalizations.

Specifically, we consider a collection of systems where is the state, is a scalar control, and is a vector of constant parameters varying on a hypercube containing a nominal parameter vector . Each system can be reduced to a scalar phase model , where the natural frequency and PRC depend on the parameter vector . In order to design a control that entrains the ensemble for all , we approximate it by , where is the nominal PRC.

Our strategy is to derive a minimum energy periodic control signal that guarantees entrainment for each system in the ensemble of oscillators