(* Setting the intial conditions *)

A = 3; 
ϵ = 0.2;
α = 3; 
γ = 0.05;   
v0 = 0.4; 
w0 = 0.4;
T = 25;
n = 50;
m = 50;

minv=-0.448775;

v0pace=0.7;
γpace=0.1;

kick = 1;

K = 1;

clustersize=Table[i,{i,1,10}];
excite = Table[z,{i,1,n},{j,1,m}];
stp = 0.25;

(* Defining the perturbed cell *)

kickleftN = Floor[(n-1)/2];
kickrightN = Floor[(n+3)/2];	
kickcellN = Floor[(n+1)/2];

kickleftM = Floor[(m-1)/2];
kickrightM = Floor[(m+3)/2];	
kickcellM = Floor[(m+1)/2];

(* Finding the intersection of the two nullclines *)

initial     = Solve[wi == A*vi*(vi-α)*(1-vi)-w0 && wi == (vi-v0)/γ, {vi,wi}, Reals];
initialline = Solve[wj == A*vj*(vj-α)*(1-vj)-w0 && wj == wi/.initial[[1]],{vj,wj},Reals];
localmax    = FindMaxValue[A*vi*(vi-α)*(1-vi)-w0, {vi,vj/.initialline[[3]]}];

startv = Evaluate[initial[[1,1,2]]];
startw = Evaluate[initial[[1,2,2]]];

(* Cycling through the Step Interval, the probability of connections and iterations *)

Monitor[Monitor[Monitor[
For[StpInt = 100,StpInt≥20,StpInt=StpInt-20,
For[prob=0,prob≤1,prob=prob+0.1,
For[JJ=1,JJ≤100,JJ++,

RC = RandomChoice[{1-prob,prob} -> {0,1},{2n+1,m+1}];

Cap=Partition[Flatten[Table[1,{i,1,n},{j,1,m}]],n];
ϵ=Partition[Flatten[Table[0.2,{i,1,n},{j,1,m}]],n];

fitzeqns =Join[

(* v equations *)

{v[1,1]'[t]==((A*v[1,1][t]*(v[1,1][t]-α)*(1-v[1,1][t])-w[1,1][t]-w0)/ϵ[[1,1]])

+RC[[1,1]]*K*(1/Cap[[1,1]])*
(0.6+0.4*(-Abs[((v[1,2][t]-v[1,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[1,2][t]-v[1,1][t])*55/(3.52278))]+StpInt)^2])*(v[1,2][t]-v[1,1][t])

+RC[[2,1]]*K*(1/Cap[[1,1]])*
(0.6+0.4*(-Abs[((v[2,1][t]-v[1,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[2,1][t]-v[1,1][t])*55/(3.52278))]+StpInt)^2])*
(v[2,1][t]-v[1,1][t])},

Flatten[Table[v[1,i]'[t]==((A*v[1,i][t]*(v[1,i][t]-α)*(1-v[1,i][t])-w[1,i][t]-w0)/ϵ[[1,i]])

+RC[[1,i-1]]*K*(1/Cap[[1,i]])*
(0.6+0.4*(-Abs[((v[1,i-1][t]-v[1,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[1,i-1][t]-v[1,i][t])*55/(3.52278))]+StpInt)^2])*
(v[1,i-1][t]-v[1,i][t])

+RC[[2,i]]*K*(1/Cap[[1,i]])*
(0.6+0.4*(-Abs[((v[2,i][t]-v[1,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[2,i][t]-v[1,i][t])*55/(3.52278))]+StpInt)^2])*
(v[2,i][t]-v[1,i][t])

+RC[[1,i]]*K*(1/Cap[[1,i]])*
(0.6+0.4*(-Abs[((v[1,i+1][t]-v[1,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[1,i+1][t]-v[1,i][t])*55/(3.52278))]+StpInt)^2])*
(v[1,i+1][t]-v[1,i][t]),
{i,2,m-1}]],

{v[1,m]'[t]==((A*v[1,m][t]*(v[1,m][t]-α)*(1-v[1,m][t])-w[1,m][t]-w0)/ϵ[[1,m]])

+RC[[1,m-1]]*K*(1/Cap[[1,m]])*
(0.6+0.4*(-Abs[((v[1,m-1][t]-v[1,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[1,m-1][t]-v[1,m][t])*55/(3.52278))]+StpInt)^2])*
(v[1,m-1][t]-v[1,m][t])

+RC[[2,m]]*K*(1/Cap[[1,m]])*
(0.6+0.4*(-Abs[((v[2,m][t]-v[1,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[2,m][t]-v[1,m][t])*55/(3.52278))]+StpInt)^2])*
(v[2,m][t]-v[1,m][t])},

Flatten[Table[v[i,1]'[t]==((A*v[i,1][t]*(v[i,1][t]-α)*(1-v[i,1][t])-w[i,1][t]-w0)/ϵ[[i,1]])

+RC[[2i-2,1]]*K*(1/Cap[[i,1]])*
(0.6+0.4*(-Abs[((v[i-1,1][t]-v[i,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i-1,1][t]-v[i,1][t])*55/(3.52278))]+StpInt)^2])*
(v[i-1,1][t]-v[i,1][t])

+RC[[2i-1,1]]*K*(1/Cap[[i,1]])*
(0.6+0.4*(-Abs[((v[i,2][t]-v[i,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i,2][t]-v[i,1][t])*55/(3.52278))]+StpInt)^2])*
(v[i,2][t]-v[i,1][t])

+RC[[2i,1]]*K*(1/Cap[[i,1]])*
(0.6+0.4*(-Abs[((v[i+1,1][t]-v[i,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i+1,1][t]-v[i,1][t])*55/(3.52278))]+StpInt)^2])*
(v[i+1,1][t]-v[i,1][t]),
{i,2,n-1}]],

Flatten[Table[v[i,j]'[t]==((A*v[i,j][t]*(v[i,j][t]-α)*(1-v[i,j][t])-w[i,j][t]-w0)/ϵ[[i,j]])

+RC[[2i-2,j]]*K*(1/Cap[[i,j]])*
(0.6+0.4*(-Abs[((v[i-1,j][t]-v[i,j][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i-1,j][t]-v[i,j][t])*55/(3.52278))]+StpInt)^2])*
(v[i-1,j][t]-v[i,j][t])

+RC[[2i-1,j-1]]*K*(1/Cap[[i,j]])*
(0.6+0.4*(-Abs[((v[i,j-1][t]-v[i,j][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i,j-1][t]-v[i,j][t])*55/(3.52278))]+StpInt)^2])*(v[i,j-1][t]-v[i,j][t])

+RC[[2i,j]]*K*(1/Cap[[i,j]])*
(0.6+0.4*(-Abs[((v[i+1,j][t]-v[i,j][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i+1,j][t]-v[i,j][t])*55/(3.52278))]+StpInt)^2])*(v[i+1,j][t]-v[i,j][t])

+RC[[2i-1,j]]*K*(1/Cap[[i,j]])*
(0.6+0.4*(-Abs[((v[i,j+1][t]-v[i,j][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i,j+1][t]-v[i,j][t])*55/(3.52278))]+StpInt)^2])*
(v[i,j+1][t]-v[i,j][t]),
{i,2,n-1},{j,2,m-1}]],

Flatten[Table[v[i,m]'[t]==((A*v[i,m][t]*(v[i,m][t]-α)*(1-v[i,m][t])-w[i,m][t]-w0)/ϵ[[i,m]])

+RC[[2i-2,m]]*K*(1/Cap[[i,m]])*
(0.6+0.4*(-Abs[((v[i-1,m][t]-v[i,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i-1,m][t]-v[i,m][t])*55/(3.52278))]+StpInt)^2])*
(v[i-1,m][t]-v[i,m][t])

+RC[[2i-1,m-1]]*K*(1/Cap[[i,m]])*
(0.6+0.4*(-Abs[((v[i,m-1][t]-v[i,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i,m-1][t]-v[i,m][t])*55/(3.52278))]+StpInt)^2])*
(v[i,m-1][t]-v[i,m][t])

+RC[[2i,m]]*K*(1/Cap[[i,m]])*
(0.6+0.4*(-Abs[((v[i+1,m][t]-v[i,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[i+1,m][t]-v[i,m][t])*55/(3.52278))]+StpInt)^2])*
(v[i+1,m][t]-v[i,m][t]),
{i,2,n-1}]],

{v[n,1]'[t]==((A*v[n,1][t]*(v[n,1][t]-α)*(1-v[n,1][t])-w[n,1][t]-w0)/ϵ[[n,1]])

+RC[[2n-1,1]]*K*(1/Cap[[n,1]])*
(0.6+0.4*(-Abs[((v[n,2][t]-v[n,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n,2][t]-v[n,1][t])*55/(3.52278))]+StpInt)^2])*
(v[n,2][t]-v[n,1][t])

+RC[[2n-2,1]]*K*(1/Cap[[n,1]])*
(0.6+0.4*(-Abs[((v[n-1,1][t]-v[n,1][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n-1,1][t]-v[n,1][t])*55/(3.52278))]+StpInt)^2])*(v[n-1,1][t]-v[n,1][t])},

Flatten[Table[v[n,i]'[t]==((A*v[n,i][t]*(v[n,i][t]-α)*(1-v[n,i][t])-w[n,i][t]-w0)/ϵ[[n,i]])

+RC[[2n-1,i-1]]*K*(1/Cap[[n,i]])*
(0.6+0.4*(-Abs[((v[n,i-1][t]-v[n,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n,i-1][t]-v[n,i][t])*55/(3.52278))]+StpInt)^2])*
(v[n,i-1][t]-v[n,i][t])

+ RC[[2n-2,i]]*K*(1/Cap[[n,i]])*
(0.6+0.4*(-Abs[((v[n-1,i][t]-v[n,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n-1,i][t]-v[n,i][t])*55/(3.52278))]+StpInt)^2])*
(v[n-1,i][t]-v[n,i][t])

+RC[[2n-1,i]]*K*(1/Cap[[n,i]])*
(0.6+0.4*(-Abs[((v[n,i+1][t]-v[n,i][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n,i+1][t]-v[n,i][t])*55/(3.52278))]+StpInt)^2])*(v[n,i+1][t]-v[n,i][t]),
{i,2,m-1}]],

{v[n,m]'[t]==((A*v[n,m][t]*(v[n,m][t]-α)*(1-v[n,m][t])-w[n,m][t]-w0)/ϵ[[n,m]])

+RC[[2n-1,m-1]]*K*(1/Cap[[n,m]])*
(0.6+0.4*(-Abs[((v[n,m-1][t]-v[n,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n,m-1][t]-v[n,m][t])*55/(3.52278))]+StpInt)^2])*
(v[n,m-1][t]-v[n,m][t])

+ RC[[2n-2,m]]*K*(1/Cap[[n,m]])*
(0.6+0.4*(-Abs[((v[n-1,m][t]-v[n,m][t])*55/(3.52278))]+StpInt)/Sqrt[1+(-Abs[((v[n-1,m][t]-v[n,m][t])*55/(3.52278))]+StpInt)^2])*
(v[n-1,m][t]-v[n,m][t])},

(* w equations *)

Flatten[Table[w[i,j]'[t]==v[i,j][t]-γ*w[i,j][t]-v0,{i,1,n},{j,1,kickcellM-1}]],
Flatten[Table[w[i,j]'[t]==v[i,j][t]-γ*w[i,j][t]-v0,{i,1,n},{j,kickcellM+1,m}]],
Flatten[Table[w[i,kickcellM]'[t]==v[i,kickcellM][t]-γ*w[i,kickcellM][t]-v0,{i,1,kickcellN-1}]],
Flatten[Table[w[i,kickcellM]'[t]==v[i,kickcellM][t]-γ*w[i,kickcellM][t]-v0,{i,kickcellN+1,n}]],
{w[kickcellN,kickcellM]'[t] ==v[kickcellN,kickcellM][t]-γ*w[kickcellN,kickcellM][t]-v0},

(* initial conditions *)

Flatten[Table[v[i,j][0]==startv,{i,1,kickleftN},{j,1,m}]],
Flatten[Table[v[i,j][0]==startv,{i,kickrightN,n},{j,1,m}]],
Flatten[Table[v[kickcellN,j][0]==startv,{j,1,kickleftM}]],
Flatten[Table[v[kickcellN,j][0]==startv,{j,kickrightM,m}]],
{v[kickcellN,kickcellM][0]==kick+startv},

Flatten[Table[w[i,j][0]==startw,{i,1,n},{j,1,m}]]

];

(* Solving the FHN equations *)

fitzsoln = NDSolve[fitzeqns,Join[Flatten[Table[v[i,j],{i,1,n},{j,1,m}]],Flatten[Table[w[i,j],{i,1,n},{j,1,m}]]],
{t,0,T},MaxSteps->Infinity];

(* Calculating the number of excited cells *)

Table[s=0;While[Evaluate[w[i,j][s]/.fitzsoln[[1]]]<localmax-0.5,
s=s+stp;If[s>T,Break[]]];
If[s>T,excite[[i,j]]=0,excite[[i,j]]=1],{i,1,n},{j,1,m}];

clustersize[[JJ]]=Total[excite,2]];

cluster=Table[clustersize[[JJ]],{JJ,1,10,1}];

Export["/Users/Rachel/Documents/MOAC/PhD/Modelling/Voltage Dependence/Step Function Model/Smooth Step Function/GJVoltDep_SF_RCS_n_"<>ToString[n]<>"_Kick_"<>ToString[kick]<>"_Width_"<>ToString[StpInt]<>"_Prob_"<>ToString[prob]<>".xls",cluster]]]

,JJ],prob],StpInt]