13
$\begingroup$

I'm interested in the periods of limit cycles of the Wilson-Cowan equations which have the form $$x'(t) = -x + S(ax(t) - by(t) +e)$$ $$y'(t) = -y + S(cx(t) - dy(t) + f)$$

where $$S(x) = 1 + \frac{\tanh(\frac{x}{2})}{2}$$

You can observe a limit cycle with the parameters seen in the following code:

    s[x_] := (1  + Tanh[x/2]/2);

    a = 10;
    b = 10;
    c = 10;
    d = -5;
    e = -0.75;
    f = -15;

    wc = {-x + s[(a*x) - (b*y) + e], -y + s[(c*x) - (d*y) + f]};
    T = 40;
    point = {0.77, 0.29};
    LimCyc = ParametricPlot[
       Evaluate[
        First[{x[t], y[t]} /. 
          NDSolve[{x'[t] == -x[t] + s[(a*x[t]) - (b*y[t]) + e], 
            y'[t] == -y[t] + s[(c*x[t]) - (d*y[t]) + f], 
            Thread[{x[0], y[0]} == point]}, {x, y}, {t, 0, T}]]], {t, 0, 
        T}, PlotStyle -> Red];

    Show[StreamPlot[wc, {x, 0, 2}, {y, 0, 2}, PlotRangePadding -> 0, 
      ImageSize -> {500, 500}], LimCyc]

Is there an easy way to numerically compute the period of a limit cycle for a given set of parameters?Limit cycle observed in Wilson-Cowan equations with parameters a = 10, b = 10, c = 10, d = -5, e = -0.75, f = -15

$\endgroup$
4
  • 1
    $\begingroup$ Could you check your parameter values? I don't get the same when I solve the model. Also, is b supposed to be d in the y'[t] equation? In general, folks prefer if you add your code to the question. And, could you also provide a reference to these equations? $\endgroup$
    – Chris K
    Aug 11, 2020 at 1:30
  • $\begingroup$ I made some edits! Sorry about that. Hopefully I have it right now. I'll try to clean up my code and add it soon. Here's a reference (eq. 7): ncbi.nlm.nih.gov/pmc/articles/PMC4733815 $\endgroup$
    – Cheyne
    Aug 11, 2020 at 1:42
  • $\begingroup$ thanks, that works now $\endgroup$
    – Chris K
    Aug 11, 2020 at 1:45
  • 3
    $\begingroup$ Check this prior MSE thread. You would want to work with res = NDSolveValue[{x'[t] == -x[t] + s[(a*x[t]) - (b*y[t]) + e], y'[t] == -y[t] + s[(c*x[t]) - (d*y[t]) + f], Thread[{x[0], y[0]} == point]}, {x[t], y[t]}, {t, 0, T}]; Plot[res, {t, 0, T}] (The plot is not actually needed, it's for illustration). $\endgroup$ Aug 11, 2020 at 14:33

3 Answers 3

13
$\begingroup$

Although it's primarily designed for ecological models, my EcoEvo package can help. First, you need to install it with

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]

Then, load the package and define your model:

<< EcoEvo`;

S[z_] := 1 + Tanh[z/2]/2;

SetModel[{
  Aux[x] -> {Equation :> -x + S[a x - b y + e]},
  Aux[y] -> {Equation :> -y + S[c x - d y + f]}
}]

Doublecheck that it matches your results:

a = 10; b = 10; c = 10; d = -5; e = -0.75; f = -15;

sol = EcoSim[{x -> 0.75, y -> 0.25}, 20];

Show[
 PlotEcoStreams[{x, 0, 2}, {y, 0, 2}],
 RuleListPlot[sol, PlotStyle -> Pink]
]

Mathematica graphics

Now use the final result from the simulation as an initial guess for FindEcoCycle:

ec = FindEcoCycle[FinalSlice[sol]];
PlotDynamics[ec]

Mathematica graphics

The period can be found as the final time of ec:

FinalTime[ec]
(* 5.27899 *)

As a bonus, you can calculate Floquet multipliers with EcoEigenvalues:

EcoEigenvalues[ec]
(* {3.6338*10^-7, -0.71155} *)

If you want to avoid the package, the idea is to warm up the simulation, look for a maximum in one variable (say x), take a tiny step further, then use WhenEvent to look for when you return to that point. There's also a Method using FindRoot.

$\endgroup$
4
  • $\begingroup$ Thanks for your help. Would the EcoEvo package also work for finding the period of a limit cycle in three dimensional model? $\endgroup$
    – Cheyne
    Aug 11, 2020 at 2:21
  • 2
    $\begingroup$ @Cheyne There's no limit on the number of equations, so maybe! $\endgroup$
    – Chris K
    Aug 11, 2020 at 2:22
  • $\begingroup$ @Chris K, I think there has been an update so , This won't work"SetModel[{ Aux[x] -> {Equation :> -x[t] + S[a x[t] - b y[t] + e]}, Aux[y] -> {Equation :> -y[t] + S[c x[t] - d y[t] + f]} }]"-- need to change to "SetModel[{ Aux[x] -> {Equation :> -x + S[a x - b y + e]}, Aux[y] -> {Equation :> -y + S[c x[t] - d y + f]} }]". Thanks for "EcoEvo" $\endgroup$ May 28, 2021 at 8:51
  • $\begingroup$ @ChrisK, Would you suggest how to implement above method if there was an explicit time dependent factor " x'[t] == -x[t] + s[(ax[t]) - (by[t]) + e] " instead of this if it were " x'[t] == -x[t](1- Sin[t] Exp[-t ]) + s[(ax[t]) - (by[t]) + e] ". It would be very kind of you , if you could show us how to implement your method for explicit time dependent differential equation . $\endgroup$ May 28, 2021 at 11:04
12
$\begingroup$

Here is a simple approach to get the period of the unknown limit cycle. The idea is to approximate the limitcycle by a circle (1st harmonic) around the mean of the limitcycle:

solution NDSolve

XY = NDSolveValue[{x'[t] == -x[t] + s[(a*x[t]) - (b*y[t]) + e],y'[t] == -y[t] + s[(c*x[t]) - (d*y[t]) + f],Thread[{x[0], y[0]} == point]}, {x, y}, {t, 0, T}]

some data of the last points

txy = Table[ { t , Norm[ Through[XY[t]]] } , {t,Subdivide[T/2, T, 100]}];

Fit the circle

{m1, m2} = NIntegrate[Through[XY[t]], {t, T/2, T}]/(T/2);
mod = NonlinearModelFit[txy, {Norm[{m1, m2} +r {Cos[2 Pi t/T1 - \[Alpha]1], Sin[2 Pi t/T1 - \[Alpha]1]}],r > 0}, { r, T1, \[Alpha]1}, t, Method -> "NMinimize"]
mod["BestFitParameters"]
(*{r -> 0.406525, T1 -> 5.28612, \[Alpha]1 -> 2.39255}*)

the period of the limitcycle T1 -> 5.28612

check result

Plot[ Evaluate[Through[XY[t]]] , {t, T/2, T},GridLines ->Evaluate[{{T - T1, T}, None} /. mod["BestFitParameters"]]]

enter image description here

$\endgroup$
2
  • $\begingroup$ Clever approach! $\endgroup$
    – Chris K
    Aug 11, 2020 at 13:16
  • 1
    $\begingroup$ @ChrisK Thanks, applying Fourier would have been much more extensive! $\endgroup$ Aug 11, 2020 at 13:22
8
$\begingroup$

I'll elaborate on the method indicated by @ChrisK involving use of WhenEvent to find a pair of maxima. Here I find a bunch of such pairs and take differences. It will be clear that they converge.

s[x_] := (1 + Tanh[x/2]/2);
{a, b, c, d, e, f} = {10, 10, 10, -5, -0.75, -15};
T = 40;
point = {0.77, 0.29};

We find both max and min values un y[t] (could also do this for x[t] but one suffices). THis is done by recording values t for which y'[t] vanishes.

extrema = 
  Reap[NDSolveValue[{x'[t] == -x[t] + s[(a*x[t]) - (b*y[t]) + e], 
      y'[t] == -y[t] + s[(c*x[t]) - (d*y[t]) + f], 
      Thread[{x[0], y[0]} == point], 
      WhenEvent[y'[t] == 0, Sow[t]]}, {x[t], y[t]}, {t, 0, 3 T}]][[2, 
    1]];

We want to go from peaks to peaks and valleys to valleys so we find time differences between extrema located two apart.

Differences[Partition[extrema, 2]]

(* Out[457]= {{5.38632, 5.29292}, {5.2813, 5.27931}, {5.27904, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}, {5.27899, 5.27899}, {5.27899, 5.27899}, {5.27899, 
  5.27899}} *)

And 5.27899 falls out as the period.

$\endgroup$
3
  • 2
    $\begingroup$ Nice simple illustration of the general idea! Just a note that while this approach works well for this example, but as a general approach it can be tricked if there are multiple local maxima within one period. In that case, you should first find a global maximum, then use WhenEvent to detect when the system returns there. $\endgroup$
    – Chris K
    Aug 13, 2020 at 16:23
  • 1
    $\begingroup$ @ChrisK Right, this method was assuming that the limit cycle is convex (weaker assumptions would suffice though). I probably should have noted that. $\endgroup$ Aug 13, 2020 at 16:46
  • $\begingroup$ It has its pitfalls, but it's nice whenever it works well. $\endgroup$ Aug 14, 2020 at 1:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.