I am trying to simulate the synaptic dynamics in a Feed-Forward system. The system is composed of a ring of $N_1$ presynaptic neurons and a single postsynaptic neuron. The presynaptic neurons are oscillating at the same frequency $\nu_1$ and are isotropically distributed across the ring. Each presynaptic neuron is described by its phase $\phi_1(i)=\frac{2 \pi i}{N_1}$.

The set of differential equations that dictate the dynamics is \begin{equation} \frac{1}{\lambda} \dot{w}(\phi_1(i), t)=\bar{w}(t)F_0(w(\phi_1(i),t))+\tilde{w}(t)F_1(w(\phi_1(i),t)). \end{equation} Where $ 1\leq i\leq N_1$, $w$ is the synaptic weight and the rest of the parameters are defined as follows \begin{equation} \bar{w}(t)=\frac{1}{N_1}\sum_{i=1}^{N_1}w(\phi_1(i),t), \end{equation} \begin{equation} \tilde{w}(t) e^{i \psi(t)}=\frac{1}{N_1}\sum_{i=1}^{N_1}w(\phi_1(i),t)e^{i \phi_1(i)}, \end{equation} \begin{equation} F_0(w(\phi_1(i),t))=D_1^2(f_+(w(\phi_1(i),t))-f_-(w(\phi_1(i),t))), \end{equation} \begin{equation} F_1(w(\phi_1(i),t))=\frac{A_1^2}{2}\tilde{K_1}\cos(\phi_1(i)-\Omega_1-\nu_1 d -\psi(t))(f_+(w(\phi_1(i),t))-f_-(w(\phi_1(i),t))). \end{equation} All the parameters that do not described as functions of time are just numbers.

I am trying to solve the first equation simultaneously for all the neurons. I did that in several ways and methods using NDSolve but no matter which method I choose I get an error. In some cases, I don't get but I know that I get a wrong result. Disintegrating the problem led me to conclude that my errors are probably due to $\psi(t)$. It seems to evaluate imaginary numbers and simply crashes or ignores them (depends on the method). I also tried to use both the function $Arg$ and $ArcTan$ to define $\psi(t)$ but it didn't help. Here is one version of my code:

(* General definitions *)

τ = 20(*ms*);
λ = 1 10^-4(*ms*);
d = 1(*ms*);
μ = 0.1;
α = 1;
N1 = 120;
D1 = 10(*Sp s^-1*) ;
A1 = 10(*Sp s^-1*);
ν1 = 2 π 6  10^-3(*ms^-1*);
fplus[w_, μ_] := (1 - w)^μ;
fminus[w_, μ_, α_] := α w^μ;
Ktilda1 := 1/Sqrt[ν1^2 τ^2 + 1];
Ω1plus = ArcTan[ν1 τ];

w1Arr[t_] := MapThread[w1, {Table[i, {i, 1, N1}], Table[t, {i, 1, N1}]}]
ϕ1 = Table[2 π i/N1, {i, N1}];

w1bar[t_] := (1/N1)Sum[w1Arr[t][[k]], {k, 1, N1}];

w1tilda[t_] := Sum[Abs[w1Arr[t][[k]]] E^(I ϕ1[[k]]), {k, 1, N1}]/N1;

w1tildaArg[t_] := Arg[w1tilda[t] // Simplify];

w1tildaAbs[t_] := Sqrt[w1tilda[t] Conjugate[w1tilda[t]]];

w1tildaArgnew[t_] := 
    Re[(1/N1)Sum[Abs[w1Arr[t][[k]]] E^(I ϕ1[[k]]), {k, 1, N1}]], 
    Im[(1/N1)Sum[Abs[w1Arr[t][[k]]] E^(I ϕ1[[k]]), {k, 1, N1}]]];

(*parametric functions*)

F01[t_, μ_, α_, i_] := 
  D1^2 (fplus[w1Arr[t][[i]], μ] - fminus[w1Arr[t][[i]], μ, α]) /. Null -> Sequence[];

  F11[t_, \[Mu]_, \[Alpha]_, i_] := 
  A1^2/2 (Ktilda1 fplus[
        w1Arr[t][[i]], \[Mu]] Cos[\[Phi]1[[i]] - \[CapitalOmega]1plus \
- \[Nu]1 d - w1tildaArgnew[t]] - 
      Ktilda1 fminus[
        w1Arr[t][[i]], \[Mu], \[Alpha]] Cos[\[Phi]1[[i]] - \
\[CapitalOmega]1plus - \[Nu]1 d - w1tildaArgnew[t]]) /. 
   Null -> Sequence[];

(* Dynamics*)

solTotal = 
       {D[w1Arr[t][[i]], t] == 
          w1bar[t]*F01[t, μ, α, i] λ + w1tildaAbs[t] F11[t, μ, α, i] λ, 
        w1Arr[0][[i]] == RandomReal[1]}, 
       {i, 1, N1}]}, 
    w1Arr[t][[1 ;; N1]], {t, 0, 10^8}, 
    Method -> {"EquationSimplification" -> "Residual"}];

Plot[Evaluate[w1Arr[t][[1 ;; 2]] /. solTotal], {t, 0, 10^4}, 
  Frame -> True, FrameLabel -> {"t", "w"}, Exclusions -> None]

In this specific case with Method -> {"EquationSimplification" -> "Residual"}, I get the wrong result. It completely ignores the right term on the RHS of the equation (the one with $\psi(t)$). By erasing the left term on the RHS it gives NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. Using other methods yield errors such as NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

Any help will be extremely appreciated.

  • $\begingroup$ When I run the code you posted it finishes with no issues or messages. It's hard to see where the problem is... $\endgroup$ – ulvi Apr 5 '18 at 5:25
  • $\begingroup$ Thank you @ulvi, I edited the code. $\endgroup$ – huo Apr 5 '18 at 13:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.