{'fraction_non_alphanumeric': 0.037935597706219674, 'fraction_numerical': 0.0171403365051358, 'mean_word_length': 3.9935494021397107, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics and optimal actuation of a three-sphere low-Reynolds-number swimmer with muscle-like arms

The three-sphere swimmer by Najafi and Golestanian is composed of three spheres connected by two arms. The case in which the swimmer can control the lengths of the two arms has been studied in detail. Here we study a variation of the model in which the swimmer's arms are constructed according to Hill's model of muscular contraction. The swimmer is able to control the tension developed in the active components of the arms. The two shape parameters and the tensions acting on the two arms are then obtained by solving a system of ordinary differential equations. We study the qualitative properties of the solutions, compute analytically their leading order approximation and compare them with numerical simulations. We also formulate and solve some optimisation problems, aimed at finding the actuation strategies maximising performance, for various performance measures. Finally, we discuss the structure of the governing equations of our microswimmers from the point of view of control theory. We show that our systems are control affine systems with drift.


Introduction
The three-sphere swimmer by Najafi and Golestanian [1,2] is a cornerstone in the literature on low Reynolds number swimming. It is composed of three spheres connected by two arms. The presence of two shape parameters (the lengths of the two arms) allows to perform periodic shape changes which are not invariant under time reversal. This is the key to beat Purcell's famous scallop theorem [3], which says that low Reynolds number swimmers cannot achieve locomotion through a reciprocal shape change. Another famous example of low Reynolds number swimmer with two shape parameters is the three-link swimmer by Purcell [3]. These two model swimmers have been studied extensively (see for example [4][5][6][7][8][9][10][11][12][13]).
In the works [1,2] the authors study a three-sphere swimmer which is able to perform periodic shape changes by controlling the lengths of the two arms. Real swimmers however are not capable of directly controlling their shape. Their shape change is generated by a complex interplay between different elements: forces generated by the swimmer through internal mechanisms, elastic properties of the swimmer's body, and forces due to interactions with the surrounding viscous fluid. The swimmer is able to control only the first of these components. For this reason it interesting to study minimal swimmers which are not capable of directly controlling their shape. An interesting work going in this direction is [14], in which the authors consider a three-sphere swimmer whose arms can generate forces thanks to a system of molecular motors and elastic elements. Another possibility is to assume that one of the two shape parameters can be controlled while the other one is driven by a passive elastic spring. This was done in [15] for the three-link swimmer and in [16] for the three-sphere swimmer.
In the present work we study a three-sphere swimmer whose arms are viscoelastic structures constructed according to Hill's active state muscle model. This model, which is very useful for computing the mechanical behaviour of muscles, includes passive elastic elements, a mechanism of viscous damping, and an active component. The swimmer is able to control the force generated by the active component. The geometry of the shape change is not known a priori and will have to be computed by solving a system of ODEs. In the first section we introduce our model and derive the ODEs governing the system. In the second section we study the qualitative properties of solutions, compute analytically the leading order term of their asymptotic expansion, present the results of numerical simulations, and study some optimmization problems. In the third section we study a variation of our model, in which one of the two muscle-like arms is replaced by a passive elastic spring.

Problem formulation
The aim of this section is to introduce the thee-sphere swimmer with musclelike arms and to obtain a system of ODEs governing its dynamics. In the first subsection we fix the notations and recall some basic facts about the threesphere swimmer. In the second subsection we assume that the forces acting through the two arms are known functions of time and shape and obtain the equations governing the evolution of the two shape parameters. In the third subsection we model the arms according to Hill's muscle model and obtain a system of ODEs governing the forces acting across each arm, the shape of the swimmer and its displacement. In the fourth subsection we write the problem in non-dimensional form. In the fifth and last subsection we consider a variation of the problem in which one of the two muscle-like arms is replaced by a passive elastic spring. Let us consider a three-sphere swimmer as in Figure 1. Let a be the radius of the spheres, and L 1 , L 2 the lengths of the two arms. We indicate with µ the dynamic viscosity of the fluid and with v i and f i respectively the velocity of sphere i and the force that this sphere exerts on the fluid, projected on the unit vector e. By using the Oseen tensor and the approximation a Li << 1 we obtain the following linear relations between forces and velocities

Three-sphere swimmer
. ( Due to Newton's third law of motion, the force exerted by the fluid on the i-th sphere is −f i . Therefore the force balance equation for the swimmer is The geometry of the system implies the following kinematic relationṡ Let us indicate with x i the position of the i-th sphere on axis corresponding to the unit vector e. We indicate with x the mean point of the three spheres, namely, The translational velocity of the swimmer is the velocity of the point x. Obviouslyẋ Using equations (1)-(4) we can show thaṫ Now we can use equations (1)-(6) to express f 1 and f 3 in terms of L 1 , L 2 ,L 1 ,L 2 . If we plug the resulting expressions into equation (9) and keep only the leading order terms in a/L j we obtaiṅ Suppose that L 1 and L 2 are periodic functions. In this case the termsL j /L j average to zero in a full swimming cycle, because they are derivatives of log(L j ). If we neglect these terms we obtain the following formulȧ which can be used instead of (10) to compute the net displacement in one period when the swimmer performs a periodic shape change. Now we consider the case of small deformations      We assume that the deformations are periodic with period T . We would like to compute the net displacement in one period to leading order in the amplitude of deformations. Let us set l := l 1 + l 2 . l is the body length of the swimmer. If we expand equation (11) and retain only the leading order terms in U i /l j we geṫ The termsU j ,U j U j ,U 1 U 2 + U 1U2 give zero when integrated between 0 and T . So the net displacement in one period l∆x is given, to leading order in U i /l j , by the following formula

TSS with assigned tensions
Let us indicate with T 1 and T 2 the tension on the tail rod and on the front rod respectively. We assume that T 1 and T 2 are known functions of time and shape, namely, T i = T i (t, L 1 , L 2 ), i = 1, 2. The equations of force balance for the three spheres are Notice that equations (15)- (17) imply that the condition of global force balance (4) is satisfied. Our aim now is to obtain two ODEs for L 1 and L 2 . From equations (15)- (17) we can easily express the forces as functions of T 1 and T 2 .
If we plug the resulting expressions into equations (1) These equations, combined with the kinematic relations (5) and (6), yield the following system of ODEs governing the evolution of L 1 and L 2

TSS with muscle-like arms
Equations (21) and (22) drive the evolution of L 1 and L 2 once T 1 and T 2 are known functions of shape and time. In this section we introduce a further element: we model the two arms of the swimmer according to Hill's active state muscle model. A schematicrepresentation of this model is shown in Figure  2. The contractile element is composed of an active component capable of generating a tension S(t) and a linear dashpot with characterisic constant B.
In addition there are two linear springs, one in parallel and one in series with the contractile element, with elastic constants k p and k s respectively. We assume that the active components of the two arms can generate tensions S 1 (t) and S 2 (t). Let us indicate with l i the sum of the rest lengths of the parallel and series springs in arm i (i = 1, 2). Tensions and lengths are related by the By rearranging the terms we finḋ By using equations (21) and (22) we can expressL 1 andL 2 as functions of T 1 , T 2 , L 1 , and L 2 . Then we can plug the resulting expressions into (10), (25), and (26). As a consequence, we can rewrite (21)-(26) and (10) as a system of five ODEs governing the evolution of x, L 1 , L 2 , T 1 , and T 2 .

Non-dimensionalization
The aim of this subsection is to write the problem in non-dimensional form. We will discover that the physical parameters influence the behaviour of the dynamics through 6 non-dimensional numbers. This approach allows to obtain more general and useful results. Let l = l 1 + l 2 be the characteristic length of the swimmer and 1/ω a characteristic time. Let us set t * = ωt, . We will indicate with the prime symbol the derivative with respect to the rescaled time t * . Using the new variables our system can be rewritten as where the dimensionless parameters are for i = 1, 2.

TSS with one muscle-like arm and one passive elastic arm
Now we replace the tail arm with a passive elastic spring. Let us indicate with l 1 the rest length and with h the elastic constant of the spring. Equations (27)-(29) and (31) remain valid. Equation (30) is replaced by where λ 1 = l 1 /l and H = h/(ωµl).

TSS with two muscle-like arms
In this section we present a detailed study of the TSS with two muscle-like arms.
For the sake of simplicity we assume that the two arms are identical. Therefore We consider the case in which the tensions developed by the active components are small, namely, with << 1 and |σ i | ≤ 1. We assume that σ 1 and σ 2 are periodic. In the first subsection we study the qualitative properties of the solutions. In the second subsection we compute the solutions analytically to leading order in .
In the third subsection we present the results of numerical simulations and in the fourth one we study some optimization problems.

Qualitative properties of the solutions
In this subsection we prove that, for small enough values of , there exists one and only one periodic orbit and that this orbit is asymptotically stable. Our analysis is based on some classical results on periodically perturbed systems which can be found in [17]. We observe that in the ODEs (27)-(31), ξ does not appear in the velocity field. This is due to the fact that the response of the system is invariant under translations. As a consequence, the differential problem for ξ is decoupled from the rest. So we can restrict our attention to the four-dimensional system for Λ 1 , Λ 2 , τ 1 , τ 2 . Instead of Λ 1 and Λ 2 we use the variables u 1 = Λ 1 − λ 1 and and (41) The unperturbed system Y = f (Y ) has the equilibrium point Y = 0. The variational system with respect to this equilibrium point (see [17], page 303) is where and Now we would like to prove that all the eigenvalues of A have a negative real part. If this is true then all the characteristic multipliers of system (42) are in modulus strictly less than one. As a consequence we can apply theorems 6.1.1 and 6.1.3 of [17] and conclude that, for small enough values of , there exists one and only one periodic orbit. Moreover, this orbit is asymptotically stable. The characteristic polynomial of A is Let us indicate with a j , j = 0, 1, 2, 3 the coefficients, so that First of all we observe that a j > 0 for every j. This follows from the fact that, if α is small enough, then Q < 0, R > 0, |Q| < R. Having all the coefficients with the same sign is a necessary condition for p A (x) to have only roots with negative real part. The assumption that α is small is a natural one, because all of our analysis is based on the hypothesis that the radius of the spheres is small compared to the lengths of the arms. The table associated to the polynomial by means of the Routh method is Let us study the signs of b 3 and c 2 . Notice that Now we notice that, since |Q| < R, we havẽ As a consequence, Notice that sign(c 2 ) = sign(a 1 b 3 − a 3 b 2 ). By using (56) we see that So b 3 and c 2 are both positive. Therefore we see from table (48) that there are three roots with negative real part. Since the determinant of A is positive, also the fourth root must have negative real part.

Asymptotic expansions
In this section we compute analytic expressions for the solutions, to leading order in . Let us assume that with |ã j | + |b j | ≤ 1, j = 1, 2. We express the solutions in power series as follows Let us set X = (u 2 ) tr . The problem satisfied by X is The solution of this problem can be computed with the help of Duhamel's formula However, since the exponential of A is difficult to compute, we choose a different strategy. We know from the previous subsection that the solution converges to a periodic orbit and we would like to study the stationary regime of the system. So we are not really interested in the solution of (64). Rather, we would like to compute the periodic orbit. This problem can be reduced to the problem of solving a linear system. First of all we notice that g can be rewritten as whereĝ s = (0, 0, J sã1 , J sã2 ) tr andĝ c = (0, 0, J sb1 , J sb2 ) tr . We look for a solution in the following form Since the equality must hold for every time t * , we obtain the following linear system We can rewrite the problem in compact form Now suppose that namely,X s = (a 1 , a 2 , c 1 , c 2 ) tr andX c = (b 1 , b 2 , d 1 , d 2 ) tr . The linear problem (70) corresponds to the following set of equations The first four equations allow us to express a 1 , a 2 , b 1 , and b 2 as linear combinations of c 1 , c 2 , d 1 , and d 2 . If we plug the resulting expressions into the last four equations we obtain where So now our problem is to check if M is invertible and to compute M −1 . Let us set Straightforward computations show that is strictly positive. So we can compute the inverse of M , and the result is where Now we can compute the coefficients The remaining coefficients can be computed using (75)-(78). At steady state the shape of the system evolves along a closed loop in the configuration space. Figure 3 shows a comparison between the loops obtained through numerical simulations and the ones corresponding to the leading order approximation.  Table 1. The different loops, from the smallest to the largest, correspond to = 0.1, = 0.3, = 0.5, and = 0.7. We see from the picture that the error becomes vanishingly small for small values of .

Numerical simulations
In this section we present the results of some numerical simulations. We used MATLAB ode45 procedure, which consists of a Runge-Kutta integration scheme with adaptive step size. The values of the dimensionless parameters used in the simulations are shown in Table 1. Table 1: Values of the parameters used in the numerical simulations.
In section 2.1 we proved that the system converges to a periodic orbit. This is confirmed by Figure 4 in which we see that, after relaxation, Λ 1 and Λ 2 evolve along a closed loop. Figure 5 shows the evolution of ξ. In this simulation (ã 1 ,b 1 ) is proportional to (ã 2 ,b 2 ) (see Table 1). This means that the tensions developed in the active components of the two arms are in phase with each other. If this sinchronization of the active components caused a sinchronization in the length change of the two arms, we would obtain a reciprocal shape change, thus no net motion in view of the scallop theorem. As we can see from figures 4 and 5, this is not the case. This property is very interesting because it implies that the swimmer can move even if the two active components are stimulated at the same frequency. We will go back to this in the following section, where this property of the system will emerge from the formula for the leading order term of the net displacement in one period. In Figure 6 we see the evolution of τ 1 and τ 2 : also in this case the behaviour at steady state is periodic, in agreement with the results of subsection 2.1.  Figure 4: TSS with two mucle-like arms: evolution of the two shape parameters Λ 1 and Λ 2 . We see that the system converges to a closed loop, in agreement with the results of section 2.1.

Optimization
In this subsection we study some optimization problems. We consider three performance measures: net displacement in one period, work per travelled distance,  Figure 5: TSS with two muscle-like arms: plot of ξ as a function of the normalized time t * /2π. We observe that in this case the net displacement in one period is negative. and Lighthill's efficiency. The value of the performance measures is determined by the dimensionless parameters that appear in our system of ODEs. These parameters depend on the swimmer's body length l, on the actuation frequency ω, on the fluid viscosity µ, and on the physical constants k s , k p , h, B characterizing the viscoelastic arms. We consider the case in which the fluid and the swimmer are given, so that the only parameter which can be varied is the actuation frequency ω. We would like to study optimization of the different performance measures with respect to ω. Our strategy is to compute the leading order approximation of the performance measures and use the resulting expressions to study our optimization problem. Then we will compare the optimality results obtained through the leading order approximation with the outcome of numerical simulations.
The first performance measure we consider is ∆x, the net displacement per period in units of body length. This quantity can be computed using formula (14). Through a simple change of integration variable we see that By using the leading order expressions computed above we find the following leading order approximation for ∆x By using equations (75)-(78) we obtain Now by using the expressions for c 1 , c 2 , d 1 and d 2 computed in subsection 2.2 we obtain (99) There are some interesting facts to point out about this formula. First of all we notice that if we change σ 1 with σ 2 and viceversa the displacement changes sign, as expected for symmetry reasons. Secondly, we observe that even in the case in which only one of the two arms is activated we obtain a non-zero displacement.
Finally, let us comment on the case in which the two active components are sinchronized. Suppose there exists a constant η > 0 such that In this case one might expect a sinchronization of the length change of the two arms, leading to a reciprocal shape change which produces no net motion.
However, our asymptotic analysis shows that is not the case provided that η = 1: from formula (99) we obtain that the leading order term of ∆x is Now we compute the mechanical work µωl 3 W done by the active components in one period. The power expenditure µω 2 l 3 P is It follows that Using the leading order approximations computed in subsection 2.2 we find that The other two performance measures we consider are are the mechanical work per travelled distance and Lighthill's efficiency Since we have already computed leading order expressions for ∆x and W , we can easily compute leading order expression for ζ and η. Now we would like to optimize the leading order terms of the different performance measures with respect to the actuation frequency ω. First of all we will show that this problem is well-posed, namely, that each performance measure admits an optimal value of ω. Secondly, we will show that the optimality results obtained by studying the leading order approximations are in good agreement with the results of numerical simulations.
So the leading order expression for ∆x vanishes when ω goes to zero and to +∞. This implies that there exists ω ∆x ∈ (0, +∞) which optimizes this performance measure. The situation is similar for the other performance measures: From these asymptotic regimes we deduce the existence of two optimal frequencies ω ζ , ω η ∈ (0, +∞). Figures 7, 8, and 9 show plots of the different performance measures as functions of ω. Notice the existence of the optimal frequencies and the very good agreement between leading order approximations and numerical simulations.   Table 1.  Table 1.

TSS with one muscle-like arm and one passive elastic arm
In this section we study the TSS with one muscle-like arm and one passive elastic arm. We assume that the active component of the muscle-like arm generates a tension with << 1 and |σ| ≤ 1. We assume that σ is a periodic functionof time.
In the first subsection we study the qualitative properties of the solutions. In the second subsection we compute the the leading order approximation of the solutions in the asymptotic regime << 1. In the third subsection we present the results of numerical simulations. In the fourth subsection we study some optimization problems.  Table 1.

Qualitative properties of the solutions
For the three-sphere swimmer with one muscle-like arm and one passive elastic arm τ 1 is given by equation (33). So we can restrict our attention to the threedimensional system of ODEs governing the evolution of the two shape variables and τ 2 . The state of the system is Y = (u 1 , u 2 , τ 2 ) tr and the problem has the form where f (Y ) is given by and (110) The unperturbed system Y = f (Y ) has the equilibrium point Y = 0. The variational system with respect to this equilibrium point (see [17], page 303) is where The characteristic polynomial of B is (116) If we prove that all the eigenvalues of B have negative real part then we can argue as in section 2.1 and conclude that, for small enough, there exists one and only one periodic orbit which is asymptotically stable. Let us call a j , j = 0, 1, 2 the coefficients of the characteristic polynomial, so that If α is small enough then P < 0, Q < 0, R > 0, |P | < R, and |Q| < R. So under this hypothesis all the coefficients of the polynomial are positive. The table associated to the polynomial by means of the Routh method is Notice that It follows that b 2 > 0. So by applying Routh's criterion we see from table (118) that there are two roots with negative real part. Since the determinant of B is negative, also the third root must have negative real part.

Asymptotic expansions
Now we would like to compute analytically the leading order term of the solutions. Let us express the solutions in power series as follows Let us set X = (u 2 ) tr . X(t * ) is the solution of The solution of this problem can be computed with the help of Duhamel's formula. Since we are interested in the periodic behaviour of the solution at steady state, we would like to compute the periodic orbit of the system. We assume thatσ (t * ) =ã sin(t * ) +b cos(t * ) .
We can write g as follows whereĝ s = (0, 0, J sã ) tr andĝ c = (0, 0, J sb ) tr . We look for the periodic solution in the following form withX s ,X c ∈ R 3 . So our differential equation becomes −X c sin(t * ) +X s cos(t * ) = (BX s +ĝ s ) sin(t * ) + (BX c +ĝ c ) cos(t * ) . (128) Since the equality must hold for every time t * , we obtain the following linear system We can rewrite the system in compact form Now suppose that namely,X s = (a 1 , a 2 , c) tr andX c = (b 1 , b 2 , d) tr . Let us rewrite (130) explicitly P Ha 1 + Rc + b 1 = 0 (134) From the first four equations we obtain If we plug these expressions into the last two equations we obtain the following problem for c and d where and The linear problem (130) has a solution if and only if (185) has a solution. So we need to check if N is invertible. We notice that hence M is invertible. The inverse of M is given by Now c and d can be computed as follows Once we know c and d we can compute a 1 , a 2 , b 1 , b 2 by using equations (140)-(143). At steady state the two shape variables evolve along a closed loop in the configuration space. Figure 10 shows a comparison between the closed curves obtained through numerical simulations and the ones corresponding to the leading order approximation.  Table 2. The different loops, from the smallest to the largest, correspond to = 0.3, = 0.5, = 0.7, and = 0.9. We can see from this picture that the error becomes vanishingly small for small values of .

Numerical simulations
In this section we present the results of some numerical simulations, obtained by choosing the dimensionless parameters as in Table 2. Figure 11 shows the J s J p K H αãb 4 3 2 5 0.1 1 0 0.7 Table 2: Values of the parameters used in the numerical simulations. evolution of the system in the shape space. We observe that, after relaxation, Λ 1 and Λ 2 evolve along a closed loop. In Figure 12 we see the evolution of ξ.
We notice that the net displacement in one period has a negative sign, which means that the object swims with the passive arm ahead. This is a general property of the TSS with passive elastic tail. Given any periodic deformation of the active arm, it can be shown that the corresponding net displacement in one period has a negative sign [16]. Figure 13 shows a plot of τ 1 and τ 2 as functions of time. The behaviour at steady state is periodic, as expected from the results of subsection 3.1.

Optimization
In this subsection we study some optimization problems for the three-sphere swimmer with one muscle-like arm and one passive elastic arm. The optimality measures we consider are the net displacement in one period ∆x, the mechanical work per travelled distance ζ and Lighthill's efficiency η. Like in section 2.4, we assume that the fluid and swimmer are given, so that the only parameter which can be varied is the actuation frequency ω. First we study our optimization problem using the leading order approximation of the different performance measures and then we compare the results with the outcome of numerical simulations.
We start by computing the leading order expressions for the different performance measures. Formula (97) for the leading order approximation of ∆x remains valid. Let us consider the mechanical work µωl 3 W 2 done by the active component of the muscle-like arm in one period. The power expenditure is

It follows that
So by using the results of section 3.2 we obtain The work per travelled distance η and Lighthill's efficiency η are defined as in section 2.4, and the corresponding leading order approximations can be easily obtained once we know the approximation of ∆x and W . Using the formulas computed in section 3.2 it is not difficult to see that the coefficients c, d, a i , b i , i = 1, 2 are O(1/ω) when ω → +∞ and O(ω) when ω → 0. Therefore, the asymptotic behaviour of the leading order terms of the different performance measures is exactly the same as the one observed in section 2.4. This implies existence of three optimal frequencies ω ∆x , ω ζ , and ω η . These theoretical predictions are confirmed by Figures 14, 15, and 16. In these figures we observe the existence of the three optimal frequencies ω ∆x , ω ζ , ω η and the high accuracy of the leading order approximations.  Table 2.

Conclusions
We studied the dynamics of the three-sphere swimmer with muscle-like arms. We assumed that the forces generated by the swimmer in the active components of the arms have intensity and vary periodically with frequency ω. We showed that the two shape parameters and the forces acting across the arms evolve according to a system of ODEs. We proved that the solutions converge to a periodic orbit. Under the assumption that << 1, we computed the leading order approximation of the solutions at steady state. Then we studied some optimization problems. We considered three different performance measures: net displacement in one period, work per travelled distance, and Lighthill's efficiency. We studied optimization of these quantities with respect to the actuation frequency ω. We computed leading order approximations of the different performance measures. Using these approximations we showed that each performance measure admits an optimal frequency. In addition we showed that the optimality results obtained through leading order approximations are in very good  Table 2. agreement with the outcome of numerical simulations. Then we introduced the three-sphere swimmer with one muscle-like arm and one passive elastic arm. We studied this model through the same type of analysis done for the three-sphere swimmer with two muscle-like arms.

A Asymptotic expansions when the forces in the active components are generic periodic functions A.1 TSS with two muscle-like arms
In section 2.2 we computed the leading order approximation of the solution in the case in whichσ 1 andσ 2 are given by equations (58) and (59). Now we would like to consider the case in whichσ 1 andσ 2 are two generic period functions  Table  2.
2 ) tr be the leading order term of the solution. X satisfies the differential equation where A and g are defined as in (43) and (41) respectively. Notice that g can be written as whereĝ s,j = (0, 0, J sã1,j , J sã2,j ) tr andĝ c,j = (0, 0, J sb1,j , J sb2,j ) tr . We are interested in computing the periodic orbit of the system, so we look for a solution in the following form X(t * ) = +∞ j=1X s,j sin(jt * ) +X c,j cos(jt * ) , withX s,j ,X c,j ∈ R 4 for every j. It is not difficult to see that the ODE (156) is equivalent to the set of linear systems for j ∈ N. Now suppose thatX We may rewrite (159) as follows Equations (162)-(165) allow us to express a 1,j , a 2,j , b 1,j , and b 2,j as linear combinations of c 1,j , c 2 j, d 1,j , and d 2,j . If we plug the resulting expressions into the last four equations we obtain the linear problems where (171) It is not difficult to check that M (j) is invertible for every j ∈ N: the computations are similar to the ones done in section 2.2 for the matrix M . So we can conclude our computation with the help of the matrices M (j) −1 , j ∈ N.

A.2 TSS with one muscle-like arm and one passive elastic arm
In this subsection we would like to generalize the results of subsection 3.2 by considering a more general periodic functionσ. We assume thatσ is a periodic function with period 2π. We expand this function in Fourier series σ(t * ) = +∞ j=1ã j sin(jt * ) +b j cos(jt * ) .
Let X = (u 1 , u 2 , τ 2 ) tr be the leading order term of the solution. X satisfies the ODE X = BX + g(t * ) , where B and g are defined by equations (112) and (110) respectively. Notice that g can be written as g(t * ) = +∞ j=1ĝ s,j sin(jt * ) +ĝ c,j cos(jt * ) , whereĝ s,j = (0, 0, J sãj ) tr andĝ c,j = (0, 0, J sbj ) tr . We would like to compute the periodic orbit of the system, so we look for a solution in the form X(t * ) = +∞ j=1X s,j sin(jt * ) +X c,j cos(jt * ) , withX s,j ,X c,j ∈ R 3 for every j. The ODE (173) is equivalent to the set of linear systems for j ∈ N. Now suppose thatX Using the first four equations we can express a 1,j , a 2,j , b 1,j , and b 2,j as linear combinations of c j and d j . Then we can plug the resulting expressions into the last to equations. The result is the following linear problem where and N 1 (j) = (J p + J s )(j 2 + P 2 H 2 ) + KH 2 (P R 2 − P 2 Q) + KHJ p R 2 /j 2 − J 2 KQ (187) N 2 (j) = j 3 + jH 2 P 2 + jKHR 2 + KH 2 J p (P 2 Q − P R 2 ) + j 2 QKJ p .
A simple computation shows that det(N (j)) > 0 for every j. So with the help of the matrices N (j) −1 we can compute the coefficients of our leading order approximation.