(* Cycling through perturbation values, probability of connection, and iterations *)

Monitor[Monitor[Monitor[
For[kick=1,kick≤2.5,kick=kick+0.5,
For[prob=0,prob≤1,prob=prob+0.1,
For[JJ=1,JJ≤100,JJ++,

n=25;
m=25;

(* Post Day-18 data *)

var = Partition[RandomVariate[NormalDistribution[(-45/-47.234706),Abs[10/-45]^2],n m],n];

(*Day-18 data*)

(*var = Partition[RandomVariate[NormalDistribution[1,0.161463468^2],n m],n];*)

(*Day-15 data*)

(*var = Partition[RandomVariate[NormalDistribution[1.04566289,0.0873686^2],n m],n];*)

(*Pre Day-15 data *)

(*var = Partition[RandomVariate[NormalDistribution[0.0955049,0.0349263^2],n m],n];*)

(* Parameter values where w0, γ are scaled according to the RMP *)

A=3;
α=3;(*Changes height of cubic nullcline peak*)
γ=Table[var[[i,j]]*0.05,{i,1,n},{j,1,m}];(*1/γ is the gradient of the straight-line nullcline*)
v0=0.4;(*-1/γ v0 is the y intercept of the straight-line nullcline*)
w0=Table[var[[i,j]]*0.4,{i,1,n},{j,1,m}];
ϵ = 0.2;
T=25;

(* Intersection of nullclines *)

initial=Table[Solve[wi==A*vi*(vi-α)*(1-vi)-w0[[i,j]]&&wi==(vi-v0[[i,j]])/γ[[i,j]],{vi,wi},Reals],{i,1,n},{j,1,m}];
initialline=Table[Solve[wj==A*vj*(vj-α)*(1-vj)-w0[[i,j]]&&wj==wi/.initial[[i,j,1,2]],{vj,wj},Reals],{i,1,n},{j,1,m}];
localmax=Table[FindMaxValue[A*vi*(vi-α)*(1-vi)-w0[[i,j]],{vi,vj/.initialline[[i,j,3]]},AccuracyGoal->4,PrecisionGoal->4],{i,1,n},{j,1,m}];

(* Defining 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];
minscale=-0.25;
maxscale=1.2;

(* Initialising connection matrix *)

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

K = 1; (* Coupling strength *)

(* Solving FHN equations *)

fitz=NDSolve[
Join[
{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*(v[1,2][t]-v[1,1][t]))
+RC[[2,1]]*(K*(v[2,1][t]-v[1,1][t]))},

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*(v[1,i-1][t]-v[1,i][t]))
+RC[[2,i]]*(K*(v[2,i][t]-v[1,i][t]))
+RC[[1,i]]*(K*(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*(v[1,m-1][t]-v[1,m][t]))
+RC[[2,m]]*(K*(v[2,m][t]-v[1,m][t]))},

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[[2 i-2,1]]*(K*(v[i-1,1][t]-v[i,1][t]))
+RC[[2 i-1,1]]*(K*(v[i,2][t]-v[i,1][t]))
+RC[[2 i,1]]*(K*(v[i+1,1][t]-v[i,1][t])),{i,2,n-1}],

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[[2 i-2,m]]*(K*(v[i-1,m][t]-v[i,m][t]))
+RC[[2 i-1,m-1]]*(K*(v[i,m-1][t]-v[i,m][t]))
+RC[[2 i,m]]*(K*(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[[2 n-1,1]]*(K*(v[n,2][t]-v[n,1][t]))
+RC[[2 n-2,1]]*(K*(v[n-1,1][t]-v[n,1][t]))},

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[[2 n-1,i-1]]*(K*(v[n,i-1][t]-v[n,i][t]))+
RC[[2 n-2,i]]*(K*(v[n-1,i][t]-v[n,i][t]))
+RC[[2 n-1,i]]*(K*(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[[2 n-1,m-1]]*(K*(v[n,m-1][t]-v[n,m][t]))
+RC[[2 n-2,m]]*(K*(v[n-1,m][t]-v[n,m][t]))},

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[[2 i-2,j]]*(K*(v[i-1,j][t]-v[i,j][t]))
+RC[[2 i-1,j-1]]*(K*(v[i,j-1][t]-v[i,j][t]))
+RC[[2 i,j]]*(K*(v[i+1,j][t]-v[i,j][t]))
+RC[[2 i-1,j]]*(K*(v[i,j+1][t]-v[i,j][t])),{i,2,n-1},{j,2,m-1}]],

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

Flatten[Table[v[i,j][0]==Evaluate[initial[[i,j,1,1,2]]],{i,1,kickleftN},{j,1,m}]],
Flatten[Table[v[i,j][0]==Evaluate[initial[[i,j,1,1,2]]],{i,kickrightN,n},{j,1,m}]],
Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[kickcellN,j,1,1,2]]],{j,1,kickleftM}]],
Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[kickcellN,j,1,1,2]]],{j,kickrightM,m}]],
{v[kickcellN,kickcellM][0]==kick+Evaluate[initial[[kickcellN,kickcellM,1,1,2]]]},
Flatten[Table[w[i,j][0]==Evaluate[initial[[i,j,1,2,2]]],{i,1,n},{j,1,m}]]],
Flatten[Join[Table[v[i,j],{i,1,n},{j,1,m}],
Table[w[i,j],{i,1,n},{j,1,m}]]],

{t,0,T},MaxSteps->100000];

(* Number of excited cells *)

CC=Table[Max[Last/@Level[Cases[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz[[1]]],{t,0,T},
PlotStyle->{{Purple,Thickness[Large]}},Background->White,PlotRange->{{-1,3.5},{-3,10}}]
,_Line,Infinity],{-2}]],{i,1,n},{j,1,m}];

clustersize[JJ]=Count[Positive[Flatten[CC-localmax]],True];
Cluster=Table[clustersize[JJ],{JJ,1,100,1}]];

Export["AsymmetryRMP_PostDay18_K_"<>ToString[K]<>"_Pert_"<>ToString[kick]<>"_Prob_"<>ToString[prob]<>".xls",Cluster]
]],JJ],prob],kick]





(*v0*)

(*vnull=Table[Solve[A*v*(v-α)*(1-v)-w-w0[[i,j]]==0,w],{i,1,n},{j,1,m}];

wnull=Table[Solve[v-γ[[i,j]]*w-v0[[i,j]]==0,w],{i,1,n},{j,1,m}];

nullclines=Table[Plot[{w/.wnull[[i,j]],w/.vnull[[i,j]]},{v,-0.5,4},Frame->True,PlotStyle->{{Red,Dashing[Large]},{Blue,Dashing[Large]}},Background->White,AxesLabel->{v,w},PlotRange->{{-1,4},{-5,10}}],{i,1,n},{j,1,m}];

initial=Table[Solve[wi==A*vi*(vi-α)*(1-vi)-w0[[i,j]]&&wi==(vi-v0[[i,j]])/γ[[i,j]],{vi,wi},Reals],{i,1,n},{j,1,m}];

trajectory=Table[ParametricPlot[Evaluate[{v[i,j][t],w[i,j][t]}/.fitz],{t,0,T},PlotStyle->{{Purple,Thickness[Large]}},Background->White,AxesLabel->{v,w},PlotRange->{{-0.5,3},{-3,3}}],{i,1,n},{j,1,m}];

Table[Show[nullclines[[i,j]],trajectory[[i,j]]],{i,1,n},{j,1,m}]*)

(*CC*)

(*Cluster*)