Interpreting the synchronisation of driven colloidal oscillators via the mean pair interaction

To explore the generic features of synchronisation facilitated by hydrodynamic interactions dominated by viscous forces, we investigate small systems of rowers; rowers are a highly simplified model for motile cilia, with each cilium approximated by a rigid sphere driven by a geometrically updated force. We introduce a new framework to analyse rowers, in which we average the pair-wise interaction of rowers by converting to the natural phase. In doing so the function describing the interaction becomes continuous, dramatically simplifying the system and allowing standard dynamical system techniques to be applied. Through inspection of phase portraits, we capture the broad features of rowers driven by power law forces and quantify their coupling strength. This approach is not limited to power law potentials, and can be applied to any monotonically increasing function. When implemented in systems with more rowers, the phase portraits show more diverse behaviour. By exploring the phase space, small systems of rowers can be designed to demonstrate specific features.


Introduction
Cilia are prevalent across eukaryotes. They are microscopic hair-like organelles that protrude from cells, which are separated into two broad classes: motile and primary. Primary cilia are present on almost all mammalian cells, and are involved in chemical and mechanical sensing [1]. Motile cilia beat periodically and are expressed by many species. In mammals they cover some of the epithelia, where they are critical for driving fluid movement along surfaces. In some single cell eukaryotes the periodic motion of their flagella is often their method of swimming [2]. In larger animals the cilia form carpets and are involved in several crucial processes: assisting in ova transport in the fallopian tubes, adjusting the flow field of cerebro-spinal fluid to target specific regions [3], and integral to mucous clearance, a critical aspect of defence against microbes and dust in the airways. For all these processes, synchronised dynamics are essential. Cilia are of such importance, that their loss of function is studied as both a primary and secondary symptom [4,5]. In primary cilia dyskinesia for example, anomalies in the cilia structure can inhibit their motion and ability to coordinate, leading to severe health complications [6]. Long term respiratory illness or trauma can also damage carpets of cilia or inhibit their ability to coordinate. Complications involving the motion of cilia have been noted in cases of cystic fibrosis [7], asthma [8], and more broadly in chronic obstructive pulmonary disease [9]. Given the range and severity of diseases involving cilia, many aspects of cilia, including their motion, have been the subject of research.
The coordination of cilia is thought to be mechanical, and there is evidence that hydrodynamic interaction plays an important role [10,11]. It has been explored using various models with differing levels of complexity. More complex models replicate the structure and motion of cilia using elastic rods or chains of beads [12][13][14][15][16]. These model cilia are then used to investigate the synchronisation between pairs of cilia or small clusters of cilia, before building up to large arrays or chains of cilia [15,[17][18][19]. The more complex models often focus on specifics like the transport speed and energy efficiency [16,[20][21][22], or the effect of variations in the cilium themselves on coordination [23]. Alternatively the cilia can be highly simplified, in which case each cilium is usually represented by a bead or sphere [24][25][26]. The advantage of this approach is any observed features are likely to be generic.
While there are several minimal models [27][28][29], this work focuses on the rower model [30,31]. In this model a cilium is approximated by a sphere that is driven by geometrically updated traps. Aspects of the pairwise interaction between rowers can be understood through an iterative map involving a simplified driving force, or through the eigenstates of the coupling tensor [32,33]. Both approaches have captured the shift from in-phase and anti-phase behaviour as the driving force is varied, and are reviewed in [34]. It is not feasible to extend the iterative map to systems of many rowers, but the eigenstate approach can be applied to larger collections of rowers. Through this lens the equilibrium behaviour of rowers in rings or chains was discussed in terms of eigenstate growth and decay between switches [33,35]. This culminated in predicting the steady state of planar arrays when varying their positions and alignment [36]. Here we provide an alternative framework to analyse the rower model, the average pair interaction. In calculating the average interaction the system become continuous, and the rowers are expressed as a dynamical nonlinear system. This ties into broader work involving Kuramoto oscillators and other general phase oscillations. The new continuous expression for the rower interaction rigorously explains the behaviour of rower pairs. When applied to larger systems, it becomes possible to determine the equilibrium solutions. Prediction of the steady state is a necessary part of designing systems with specific behaviour, which might inform research into engineering functional collective filament dynamics in artificial systems.

Rower model
Simple models for cilia approximate both the fluid dynamics, and the filament internal degrees of freedom, by replacing the filament with a sphere. In the rower model the bead oscillates along a direction, say x, driven by potential traps that are updated geometrically. The shape of the potential contains in a coarse-grained fashion the degrees of freedom of the complex cilium shapes and activity, and the far-field fluid dynamics can be matched to a given biological system [34]. The potential in the model can be any monotonically increasing function, but we predominantly focus on a simple power law with a kx ; r k is the trap strength and x r the distance relative to the vertex of the trap. Details of the model are represented in figure 1, with aspects of the power law potentials covered in 1(b). The position update behind the oscillations is illustrated in figure 1(a). Once a rower passes a certain threshold x s the trap is reflected and the bead reverses direction, with the threshold measured relative to the vertex. The reflection axis is chosen to create oscillations with amplitude A. The period of the oscillation is set by the driving potential and the drag in the fluid, whereas the phase of the oscillation is free and subject to drift in the presence of Brownian noise.
The one alternative to the power law potential we consider is a monotonic region of a cubic potential. Specifically we chose a trapping potential that included a non-stationary point of inflection. This ensures the potential has three stages, slow-fast-slow, rather than switching from fast-slow or slow-fast as in the power law case. The form of the potential is = - , where x R is the distance between the local minimum and local maximum, see figure 1(d). The switch point x s controls the proportion of the trapping potential that has positive curvature. The effect of increasing x s is illustrated in figure 1(e). For small x s the potential more frequently has positive curvature, marked in green. When x s is increased the potential gains sections with negative curvature, indicated by purple. When = - the driving potential is centred about the point of inflection, and has equal parts positive and negative curvature. This potential is more similar to that of Chlamydomonas when its flagella are approximated by spheres moving on the filament's centre of drag [32,37]. The trap centres are spaced d apart, so the distance between a pair of rowers is on average d, see figure 1(c). The y position of the rowers is maintained by a harmonic trap centred on 0, restricting oscillations to onedimension along the x-axis.
The rowers are coupled by hydrodynamic forces, because bead movement sets up a velocity field. The equation governing a given rower is, where f j is a stochastic noise term, F j the trapping force, and H ij represents the hydrodynamic interaction between the N rowers. From Brownian Dynamics the noise has zero mean and [38]; k B is Boltzmann's constant and T the temperature. In the simplest conditions (far field approximation, and absence of other boundaries) the hydrodynamic interaction is implemented through the Oseen tensor, The drag is g ph , where η is the viscosity of the fluid, a the radius of the rower sphere, I the 3×3 identity, andr ij and r ij are the direction and distance between rowers.

Average pair interaction
3.1. Natural phase gauge To calculate the average interaction between rowers we first transform into the natural phase frame. This transformation is often used to link Josephson junctions to Kuramoto phase oscillators, but it is not restricted to such cases [28,39,40]. The natural phase is a frame where the velocity is constant for an uncoupled rower. The natural phase accounts for any variations in velocity due to the driving force, but not for velocity changes stemming from coupling between oscillators. The transformation also reduces the parameters describing the place a rower is in its cycle to just the phase f, rather than its position and the trap direction.
The relationship between phase and position can be determined using the velocity in terms of x as determined by equation (1). We approximate by dropping the noise and only considering the direction of oscillation x. Equations (1) and (2) then simplify to: Figure 1. Rowers are beads driven by a geometrically updated force, they were designed as a highly simplified model for motile cilia and have been studied both theoretically and experimentally. A rower is driven by a potential that updates at a given threshold x s , see (a). At x s the potential is reflected and the rower changes direction; in (a) this is the switch from the dark green to light green curve. The amplitude of the resulting oscillations is A. (b) The driving potential can be any monotonically increasing function, but we focus predominantly on the simple case where the potential is a power law µ a x r , with x r the distance relative to the trap edge. (d) The only extension to power law potentials we consider is a cubic potential with a mixture of positive (marked in green) and negative (purple) curvature, i.e. near a non-stationary point of inflection. The switch point x s is measured relative to the local minimum, and x R is the distance between the minimum and maximum.
The hydrodynamic coupling in the x direction is H ij x and the trap force along this direction F j (x). To simplify future calculations, we approximate the distance between rowers (r ij ) by the distance between their trap centres ( ( ) r ij t ). The hydrodynamic coupling for a given pair is constant under this approximation, which we label μ ij , x ij t is used to indicate the trap separation in the x direction. This approximation assumes the oscillation The reparametrisation of position to the natural phase transforms a rower into a moving reference frame. For an uncoupled rower in this frame the phase velocity is constant. To determine the transformation to this frame we consider a single rower, and relate its phase to position using the chain rule, Ω is the constant phase velocity for an isolated rower, which we defined using τ the time between successive trap updates of a single rower i.e. the semi-period. We choose the natural phase to be strictly increasing. It is easier to relate the phase to position using the position relative to the trap vertex = - which is strictly decreasing with time. Given a specific trap force, the transformation between the natural phase and the position coordinates can be calculated by integrating. This transformation between phase and position is used to express the driving force in terms of phase, and so include interactions between rowers in the phase velocity. This is expanded upon in the following section.

Cycle average
When including interactions between rowers in the evolution of the natural phase, the interaction can be approximated by its average effect over a period. This approach has two key benefits: the interaction is explicitly dependent on the differences between rowers, and the interaction becomes continuous. This allows the system to be analysed in terms of fixed points and stability using standard techniques. This is a far more tractable problem than previous approaches [32,33].
The approximation relies on weak coupling between the rowers, i.e. μ ij is small, which is true in the far field we are considering. The full derivation of the phase evolution is laid out in appendix A, with the final expression, is defined as the mean of the interaction function in the natural phase frame, where the variable ϑ is a dummy used for the averaging. This expression is generic and applies to any potential, assuming the phase/ position inverse exists i.e. the potential is monotonically increasing.
When considering the evolution of the phase difference for a single pair of rowers, the even component of the interaction is the same regardless if the rower is leading or trailing ( where μ 12 =μ 21 . Calculating the effective interaction by transforming to the natural phase gauge and taking the average over a period uses equivalent steps to those taken in [28]. It is remarkable that in doing so the expression used for the interaction is no longer discontinuous. This converts the complex configuration and historydependent interactions of a rower pair synchronising to a simple one-dimensional nonlinear system. The fixed points of the system and their stability can be investigated from this perspective, where before it was interpreted in terms of eigenstates of the hydrodynamic tensor and their decay rate [33,35,36].

Representing the interaction graphically for power law potentials
The mean interaction for power law potentials can be calculated analytically, details are in appendix B.1. The resulting expression y [ ] G int depends on two variables p and C 0 , a a The parameter p depends exclusively on the trap shape α.
0, 1 to α>1 and for positive curvature. α=2 is the special case where p diverges. The divergence stems from the calculation of phase-position transformation which involves that assumed a ¹ 2. Considering the case α=2 separately, the integral still has a solution and it is still possible to solve for y [ ] G int . C 0 ensures the switch point of the trap corresponds to f=nπ. It depends on the amplitude and switch point, so in some way quantifies the variation of the force within the trap. A higher value corresponds to a potential with more similar forces and hence weaker coordination.
Examples of the interaction are shown in figure 2. In figure 2(a) p=2/7 (α=0.6), and the curvature is negative. The interaction shows that when a rower lags behind another (ψ<0), it experiences a force greater than the rower ahead. This is confirmed in 2(d), where the odd component of the interaction is plotted. The discrepancy in the force leads to the lagging rower catching the leading rower, and so in-phase motion ψ=0 is stable. Figure 2(b) is the special case where p=0 and the trap potential is linear. In this case the average interaction is an even function, emphasised in figure 2(e), and rowers experience the same force regardless if they lead or trail. Consequently there is no impetus to synchronise the rowers, as explored experimentally in [32]. In 2(c) p=−3/7 (α=1.3) and the trap has positive curvature. The average interaction in this case shows the leading rower experiencing greater force than the lagging rower. This exacerbates the separation between the rowers and pushes the system towards anti-phase. Figure 2(f) confirms the stability of the anti-phase point. In calculating the average interaction we have accounted for the known shift between anti-phase and in-phase behaviour when varying the shape of a power law potential.
Next we go further and quantify the coupling strength associated with a given potential.

Fourier series approximation of the relaxation time
To quantify the coupling strength of a given potential we use the relaxation time. Capturing the synchronisation strength increases the utility of the mean interaction approach. At this stage it is useful to express the mean interaction in terms of its Fourier series. This is a simple way to separate the odd and even contributions, and accelerates any calculations. The series will converge, because the averaging process ensured the interaction function is continuous. A side benefit of this is that it ties back to more generic work with Kuramoto oscillators and other general oscillators [39,41,42]. The sine coefficients s k are calculated using the dimensionless form of G int , and so the evolution of the phase difference is expressed as, The even sine terms s 2k were zero across all the cases we investigated. Any driving potential with a moderate s 2 term would be of interest, because we believe it is the most likely case to exhibit phase-locking where y p ¹ n . We find that the first few non-zero terms are enough to capture the details of the system.
To calculate the relaxation time we assume that the system is near equilibrium. We then linearise the system, which leads to an exponential solution. Assuming the in-phase point is stable then, sin sin 3 ...
Including more terms increases the accuracy close to the fixed point but also limits the region it is applicable, due to restrictions when linearising. It was assumed that y y » ( ) k k sin during linearisation, which completely breaks down when y p > k 2. For this reason keeping to a minimal number of terms is often more robust than including many additional terms. The first order relaxation time of the power law potentials is displayed in figure 3(a): it diverges when  p 0 reflecting the increasing difficulty to synchronise as the driving potential becomes more linear. The relaxation time also becomes longer as C 0 is increased, reflecting that for large x s the forces don't vary particularly within the trap.

Power law potentials
We verified the analytical relaxation time by comparing with the relaxation time in simulation. The pair of rowers from section 4 are simulated by updating their position using the Langevin equation in (1). The noise is implemented using the Ermack McCammon method [38]. The dimensionless noise ξ and velocity parameterṼ are, The average force over the duration of the trap á ñ 2 ) and the trap randomly oriented. Once the pair reaches equilibrium the relaxation time is calculated using the autocorrelation of the phase difference. We measure the relaxation time by fitting an exponential to the autocorrelation and recording the rate of decay. This simulation run is repeated 20 times with different initial conditions. The extreme outliers are then removed before calculating the mean and standard deviation of t r ; extreme outliers For small distances, adding terms improves the in-phase approximation but the anti-phase fit deteriorates. This discrepancy disappears at large d, and is likely due to a breakdown of the fixed distance assumption. The effect of noise is shown: high noise in grey (ξ=1.93×10 −3 ) and low noise in black (ξ=9.70×10 −4 ). The higher noise case is usually better captured by lower order approximations, because increasing the order also decreases the region where they are valid. The insets zoom in on the analytical results, with the convention that ( ) t r 1 , ( ) t r 3 , ( ) t r 5 , and ( ) t r 7 correspond to solid, dashed, dotteddashed, and dotted lines. For the in-phase case ( ) t r 13 is also included, the line with large dots, to coincide with the strongly coupled cases at low d/A. (c)The predicted relaxation time is compared with the simulated results when the trap shape is changed. The inset emphasises the decreasing effect of adding correction term, using the same convention as in (b) but only considering up to ( ) t r 5 . We set d/A=10 when changing the trap shape, to justify approximating the hydrodynamic coupling as a constant. (d)The predicted relaxation time is accurate also when considering the more complex cubic potential; here the change in trap shape is captured by the average curvature á ñ c . (e)The change stability is predicted by the mean interaction, with the analytical fixed point (orange curve) coinciding with the average phase difference y á ñfrom simulations (grey squares).
occasionally occur in weakly coupled cases, and are defined as being far outside the interquartile range, 1 with Q 1 and Q 3 the first and third quartiles. The comparison between the simulated and analytical Fourier series relaxation times are shown in figures 3(b), (c). The trap shape is varied in (c), and the simulated results closely follow the Fourier series approximation. The inset focuses on changes as the number of Fourier terms in the linearisation is increased. In both the plot and the inset the first order term is the solid line, the third is the dashed, and the fifth the dotteddashed. The simulation measurements for relaxation time are marked with the grey error bars. The bars mark the confidence interval  S t r is the measured standard deviation, and the number of trials is 20, n=20. The trap separation is set to d/A=10, to ensure the hydrodynamic coupling is approximately constant throughout the cycle, i.e. m » H x 12 12 . The relationship between rower separation and relaxation time is explored in more detail in figure 3(b). The top panel (green lines) with p=0.17 has the pair moving in-phase, and the rowers are expected to have a constant separation throughout a period. The grey and black markers represent different noise levels, grey-ξ=1.93×10 −3 , black-ξ=9.70×10 −4 . Higher noise results are often better represented by t r with fewer terms. Including more terms restricts the region about the fixed point in which the measure is valid. The restriction results from the linearisation in section 4. In the bottom panel (purple lines) p=−0.25 and the rowers are synchronising in anti-phase. When moving anti-phase the distance between rowers varies through the cycle, and for small separation the approximation m = H x 12 12 is poor. In this case we find the simulation relaxation time is higher than would be expected at small separation. The mismatch between simulation and analytical relaxation time is exacerbated by including additional terms. At larger separation the results are more akin to the in-phase case, with lower noise cases better represented by t r with more terms than the higher noise.

Chlamydomonas style cubic potential
The mean interaction can be calculated for any monotonic potential. We also consider a cubic potential that has similarities to the resistive force approximation of Chlamydomonas flagella. For consistency all the parameters of the cubic potential, including the average force, are the same as in the power law case. The noise is set to an intermediate value of the two previously used, ξ=1.45×10 −3 . In this case the trap shape is varied through the switch point x s , not α. To ensure a mixture of curvature in each trap we set

The small offset V =
A 0.03 is included to prevent the rower stagnating at the extremes, i.e. where the driving force is zero.
The analytical expression for y [ ] G int in the cubic case also exists and can be found in appendix B.2, but it and the associated parameters lack a natural interpretation. Instead we will consider the average curvature of the potential á ñ c , ò á ñ = - The relaxation time is plotted against the average curvature in figure 3(d). The theoretical values from G int [ψ] are shown by the orange lines, with the same style convention as in figure 3(c). The results from simulation are marked with grey squares. The agreement between the simulation and Fourier series relaxation time is as high for the more complex, cubic case as it is for the power law.
The fixed point as predicted by the average interaction also coincides with the equilibrium behaviour seen in simulation. The average phase difference (grey squares) and the phase difference of the stable fixed point (orange line) are plotted against á ñ c in figure 3(e). The two are in agreement for all á ñ c . At á ñ = c 0 the average interaction has no stable fixed point and the simulated rowers do not synchronise, reflected in the large error bars.
Previously the equilibrium behaviour for a potential derived from the Chlamydomonas swimming stroke could only be corroborated using simulation, and the synchronisation strength inferred through sensitivity to noise. Now with the average interaction we can predict the behaviour and the relaxation time of the potential. The relaxation time quantifies the synchronisation strength of different potentials, and can be used to compare driving forces. In this way it is possible to probe the sensitivity of a potential to small changes in shape.

The mean interaction as a predictive tool
The new approach of calculating an average interaction simplifies the interaction between a pair of rowers to a one-dimensional nonlinear equation. Including additional rowers increases the dimension of the system of equations. If restricted to a particular driving force, i.e. a particular Fourier series, we can explore the effect of varying the coupling between different rowers. In this way we can search for desired or interesting behaviour and then design the rower layout to achieve the necessary coupling. This extends the use of the average interaction beyond just understanding the synchronisation between rower pairs to a predictive tool for designing systems of rowers.
We focus here on the three rower case, which has three pair combinations and so three different coupling strengths. To reduce this, we assumed the coupling was symmetric between two of the pairs. An isosceles triangle layout is a simple geometry that would lead to this sort of dynamic, represented in figure 1(f). In this system the dynamics depend on a single parameter R μ , the ratio between the unique pair coupling μ LR and the symmetric coupling μ LC =μ RC . We can tailor the system for any value by varying d v , the vertical distance between the central rower (C) and the left and right rowers (L, R), demonstrated in figure 4(a). We used the Fourier series derived from a power law potential with p=−0.25 and C 0 =0.21.
The evolution of the phase difference between the left-right pair, y LR , and left-centre pair, ψ LC , can be calculated by extending equation (8) to three rowers, i.e Î { } i j , 1, 2, 3 . Using these equations the phase portraits of the system for a given R μ can be plotted, examples of which are shown in figure 1(g). The nullclines are marked by purple and green curves, with the fixed points indicated by orange points; unstable are empty and stable points are filled circles. If a stable limit cycle exists it is also marked in orange, but is a curve. The portraits show the saddle-node bifurcation that occurs when the coupling with rower C decreases. The phase space is toroidal and must be bounded, so a stable limit cycle forms as the fixed points 'mask' and then annihilate one another. Examples of the two stable behaviours are shown using Poincaré sections in figure 4(b). ψ LR is plotted at each point rower C completes a cycle, i.e. phase of rower C f ( c) =2π. The same parameters as in section 5.1 are used, but the noise is decreased to emphasise the limit cycle behaviour, ξ=4.85×10 −4 .
The rowers are phase-locked at the fixed points, the predicted shift in phase difference is compared to simulation in figure 4(c). The curves are the predicted values and the markers the results of simulation. The colour indicates the pair being considered, green for ψ LR and purple for ψ LC . The marker changes from a square to triangle when the simulations are displaying a limit cycle. The centre of the limit cycle is calculated by finding the midpoint of cycle and averaging. It is meaningless to average ψ LC for the limit cycle because it continually grows, see the inset. At R μ =2.4 both the limit cycle behaviour and the phase-locked state were observed. The two behaviours were considered separately and both results shown, hence the presence of the square and triangle markers. This aligns with the predictions from the phase portraits, which show a region of coexistence between the phase-locked state y á ñ » 2 LR and the stable limit cycle centred at y p á ñ = LR . The phase portraits derived using the average interaction are a highly effective predictive tool for rower systems. Here we focused on a given driving force for a three rower system, where one rower was designed to have equal coupling with the other two. The phase portraits predicted a shift from a fixed point with a set phase difference to a limit cycle. Of particular interest was the overlap of the regions where each state is stable, revealing a region of bistability. This approach is not limited to the case discussed here, but can be applied to any system with an arbitrary number of rowers. It is only limited by the ability to visualise and understand the phase space. Poincaré sections of simulations reveal a limit cycle centred on π for R μ =3.1, and a stable fixed point near ψ LR =2 for R μ =1.1. (c) This is consistent with the predictions from the phase portraits, which initially showed a stable fixed point, but as R μ is increases disappears and is replaced by a stable limit cycle. The simulation results for the fixed points are shown by square markers, with green and purple distinguishing between ψ LR and ψ LC . For the limit cycle the mean of ψ LC is irrelevant, because ψ LC continues to grow in time (see inset). The mean of the cycle for ψ LR is marked with green triangles, where the error bars are calculated using the standard deviation of the each cycles midpoint. Near R μ =2.2 both the limit cycle and phaselocked state were observed, when plotted each behaviour was considered separately. The analytical results for the fixed point are marked by the green and purple lines.

Conclusion
We have shown the mean interaction between a pair of geometrically updated oscillators captures their synchronisation properties. The structure of the average interaction's odd component accounts for the switch between in-phase and anti-phase behaviour when the trap shape is changed. This explanation involving a onedimensional nonlinear system is more rigorous than other attempts. Previously the change from in-phase to anti-phase behaviour was interpreted as a preference for the decay or growth of an eigenstate between trap updates. While this does capture the change, there is no obvious method to interpret different levels of the decay rate in terms of synchronisation strength. In contrast this new approach captures the change in the steady state, and measures the coupling strength by predicting the relaxation time. This allows direct comparison of the synchronisation strength of different potentials. Further, it ties rower systems to the larger body of work involving bifurcations and fixed point analysis, which was previously inaccessible The mean interaction technique can be applied to any monotonically increasing driving potential. Here we included a cubic potential with varying levels of positive and negative curvature. The mixture of slow and fast regions is more similar to the potentials based on resistive force theory applied to actual flagella. Through the mean interaction it is again possible to predict the equilibrium behaviour and relaxation time of these more complex potentials. As this analysis can be applied to any potential, it is possible to explore the sensitivity of the steady state and relaxation time to changes in the interaction.
The benefits of the mean interaction framework become even more apparent when applied to larger systems of rowers. With more rowers, more diverse behaviour can be exhibited. In this case the trap potential is no longer the only defining feature of the equilibrium state. The coupling relative to other rowers, the direct consequence of their physical positions, also has a key role in determining the type of behaviour. The mean interaction can be used to explore both the relative coupling as well as the effect of the trapping potential. The separation of the hydrodynamic coupling and the trap means the average interaction can be used as a predictive technique. The relative coupling can be varied to find specific behaviour, and then the rower configuration can be designed to produce the desired outcome. This can be of assistance when building blocks of rowers that are sensitive to changes of a single, control rower.
The natural phase has a constant velocity when μ ij =0, but this is no longer true when the rowers interact. The counter N s enumerates the number of trap updates that have taken place, which distinguishes between cases with x t >0 and x t <0.  f 0 shifts the phase such that f is a multiple of π when the rower reverses direction. N s is a counter that follows the same convention as before, with the switch point x s at nπ and f p p Î [ ) , 2 for x t <0. In the cubic case, the force defined in terms of the natural phase is,