ClearAll["Global`*"];
Remove["Global`*"];

Monitor[Monitor[Monitor[

(* Incrementing through perturbations, probability of connections, iterations *)

For[kick=0.5,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;


(*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];

(* Parameter values *)

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=Table[0.4,{i,1,n},{j,1,m}];(*-1/γ v0 is the y intercept of the straight-line nullcline*)
w0=Table[var[[i,j]]*0.4,{i,1,n},{j,1,m}];
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 structure *)

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

(* Day 18 *)

(*Cap = Partition[RandomVariate[LogNormalDistribution[−0.035506,0.266481^2],n m],n];
ϵ = Cap/(Exp[−0.035506-0.266481^2]);*)


(*Day-15*)

		Cap = Partition[RandomVariate[LogNormalDistribution[−0.21139,0.307329^2],n m],n];
ϵ = Cap/(Exp[−0.21139-0.307329^2]);

(* Resistance defined by capacitance and a sf *)

sf = 1;
R =sf/Cap[[1,1]];

(* 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]])/ϵ[[1,1]])
+RC[[1,1]]*((1/(R*Cap[[1,1]]))*(v[1,2][t]-v[1,1][t]))
+RC[[2,1]]*((1/(R*Cap[[1,1]]))*(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]])/ϵ[[1,i]])
+RC[[1,i-1]]*((1/(R*Cap[[1,i]]))*(v[1,i-1][t]-v[1,i][t]))
+RC[[2,i]]*((1/(R*Cap[[1,i]]))*(v[2,i][t]-v[1,i][t]))
+RC[[1,i]]*((1/(R*Cap[[1,i]]))*(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]])/ϵ[[1,m]])
+RC[[1,m-1]]*((1/(R*Cap[[1,m]]))*(v[1,m-1][t]-v[1,m][t]))
+RC[[2,m]]*((1/(R*Cap[[1,m]]))*(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]])/ϵ[[i,1]])
+RC[[2 i-2,1]]*((1/(R*Cap[[i,1]]))*(v[i-1,1][t]-v[i,1][t]))
+RC[[2 i-1,1]]*((1/(R*Cap[[i,1]]))*(v[i,2][t]-v[i,1][t]))
+RC[[2 i,1]]*((1/(R*Cap[[i,1]]))*(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]])/ϵ[[i,m]])
+RC[[2 i-2,m]]*((1/(R*Cap[[i,m]]))*(v[i-1,m][t]-v[i,m][t]))
+RC[[2 i-1,m-1]]*((1/(R*Cap[[i,m]]))*(v[i,m-1][t]-v[i,m][t]))
+RC[[2 i,m]]*((1/(R*Cap[[i,m]]))*(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]])/ϵ[[n,1]])
+RC[[2 n-1,1]]*((1/(R*Cap[[n,1]]))*(v[n,2][t]-v[n,1][t]))
+RC[[2 n-2,1]]*((1/(R*Cap[[n,1]]))*(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]])/ϵ[[n,i]])
+RC[[2 n-1,i-1]]*((1/(R*Cap[[n,i]]))*(v[n,i-1][t]-v[n,i][t]))+
RC[[2 n-2,i]]*((1/(R*Cap[[n,i]]))*(v[n-1,i][t]-v[n,i][t]))
+RC[[2 n-1,i]]*((1/(R*Cap[[n,i]]))*(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]])/ϵ[[n,m]])
+RC[[2 n-1,m-1]]*((1/(R*Cap[[n,m]]))*(v[n,m-1][t]-v[n,m][t]))
+RC[[2 n-2,m]]*((1/(R*Cap[[n,m]]))*(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]])/ϵ[[i,j]])
+RC[[2 i-2,j]]*((1/(R*Cap[[i,j]]))*(v[i-1,j][t]-v[i,j][t]))
+RC[[2 i-1,j-1]]*((1/(R*Cap[[i,j]]))*(v[i,j-1][t]-v[i,j][t]))
+RC[[2 i,j]]*((1/(R*Cap[[i,j]]))*(v[i+1,j][t]-v[i,j][t]))
+RC[[2 i-1,j]]*((1/(R*Cap[[i,j]]))*(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];

(* Calculating 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_Day15_Cap15_ScaleFactor_"<>ToString[sf]<>"_Pert_"<>ToString[kick]<>"_Prob_"<>ToString[prob]<>".xls",Cluster]
]],JJ],prob],kick]