Monitor[Monitor[Monitor[

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

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++,

(* setting parameter values *)

A = 3; 
ϵ = 0.2;
α = 3; (* Changes height of cubic nullcline peak *)
γ = 0.05;   (* 1/γ is the gradient of the straight-line nullcline *)
v0 = 0.4; (* -1/γ v0 is the y intercept of the straight-line nullcline *)
w0 = 0.4;
T = 25;
n = 25;
m = 25;


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]]}]-wi/.initial[[1]];

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

maxw1 = 0;

minscale=-0.25;
maxscale=1.2;

(* Setting the uniform distributions for connection strengths *)

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];
RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0];


(* Excluding cases where the perturbed cell is isolated *)

While[Mean[{RRh[[kickcellN,kickcellM]],RRh[[kickcellN,kickcellM+1]],RRv[[kickcellN,kickcellM]],RRv[[kickcellN+1,kickcellM]]}]==0,

RRh = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n},{m+1}];
RRv = Table[RandomChoice[{1-prob,prob}->{0,RandomReal[{0.5,10}]}],{n+1},{m}];
RRh = ReplacePart[RRh,{_,1|m+1}->0];
RRv = ReplacePart[RRv,{1|n+1,_}->0]];


(* Closed Boundary Conditions *)
(* Solving the 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)/ϵ)
+RRh[[1,2]]*(v[1,2][t]-v[1,1][t])
+RRv[[2,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)/ϵ)
+RRh[[1,i]]*(v[1,i-1][t]-v[1,i][t])
+RRv[[2,i]]*(v[2,i][t]-v[1,i][t])
+RRh[[1,i+1]]*(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)/ϵ)
+RRh[[1,m]]*(v[1,m-1][t]-v[1,m][t])
+RRv[[2,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)/ϵ)
+RRv[[i,1]]*(v[i-1,1][t]-v[i,1][t])
+RRh[[i,2]]*(v[i,2][t]-v[i,1][t])
+RRv[[i+1,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)/ϵ)
+RRv[[i,m]]*(v[i-1,m][t]-v[i,m][t])
+RRh[[i,m]]*(v[i,m-1][t]-v[i,m][t])
+RRv[[i+1,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)/ϵ)
+RRh[[n,2]]*(v[n,2][t]-v[n,1][t])
+RRv[[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)/ϵ)
+RRh[[n,i]]*(v[n,i-1][t]-v[n,i][t])
+RRv[[n,i]]*(v[n-1,i][t]-v[n,i][t])
+RRh[[n,i+1]]*(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)/ϵ)
+RRh[[n,m]]*(v[n,m-1][t]-v[n,m][t])
+RRv[[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)/ϵ)
+RRv[[i,j]]*(v[i-1,j][t]-v[i,j][t])
+RRh[[i,j]]*(v[i,j-1][t]-v[i,j][t])
+RRv[[i+1,j]]*(v[i+1,j][t]-v[i,j][t])
+RRh[[i,j+1]]*(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]-γ*w[i,j][t]-v0,{i,1,n},{j,1,m}]],

			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,1,kickleftN},{j,1,m}]],
			Flatten[Table[v[i,j][0]==Evaluate[initial[[1,1,2]]],{i,kickrightN,n},{j,1,m}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,1,kickleftM}]],
			Flatten[Table[v[kickcellN,j][0]==Evaluate[initial[[1,1,2]]],{j,kickrightM,m}]],
			{v[kickcellN,kickcellM][0]==kick+Evaluate[initial[[1,1,2]]]},
			Flatten[Table[w[i,j][0]==Evaluate[initial[[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];


(* Calculate the number of excited cells *)

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

Export["Network Topology/Isotropic_RCRS_Prob_"<>ToString[prob]<>"_Pert_"<>ToString[kick]<>".xls",Cluster]]]


,JJ],kick],prob]