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

Monitor[Monitor[

(* Cycling through connection probability and iteriations *)

For[prob=0,prob≤1,prob=prob+0.1,
For[JJ=1,JJ≤100,JJ++,

 (* parameters *)
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;
n = 25;
m = 25;

T = 200;

dt = 0.05;
tmin = 0;
tmax = T;

γpace = 0.09;
v0pace = 0.7;

K = 1; (* Coupling strength *)

(* Intersection of non-pacemaker 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]]}];

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

maxw1 = 0;

minscale=-0.4;
maxscale=3.2;

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

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

Flatten[Table[v[i,kickcellM]'[t]== ((A*v[i,kickcellM][t]*(v[i,kickcellM][t]-α)*(1-v[i,kickcellM][t])-w[i,kickcellM][t]-w0)/ϵ)+RC[[2i-2,kickcellM]]*(K*(v[i-1,kickcellM][t]-v[i,kickcellM][t]))+RC[[2i-1,kickcellM-1]]*(K*(v[i,kickcellM-1][t]-v[i,kickcellM][t]))+RC[[2i,kickcellM]]*(K*(v[i+1,kickcellM][t]-v[i,kickcellM][t]))+RC[[2i-1,kickcellM]]*(K*(v[i,kickcellM+1][t]-v[i,kickcellM][t])),{i,2,kickcellN-1}]],

Flatten[Table[v[i,kickcellM]'[t]== ((A*v[i,kickcellM][t]*(v[i,kickcellM][t]-α)*(1-v[i,kickcellM][t])-w[i,kickcellM][t]-w0)/ϵ)+RC[[2i-2,kickcellM]]*(K*(v[i-1,kickcellM][t]-v[i,kickcellM][t]))+RC[[2i-1,kickcellM-1]]*(K*(v[i,kickcellM-1][t]-v[i,kickcellM][t]))+RC[[2i,kickcellM]]*(K*(v[i+1,kickcellM][t]-v[i,kickcellM][t]))+RC[[2i-1,kickcellM]]*(K*(v[i,kickcellM+1][t]-v[i,kickcellM][t])),{i,kickcellN+1,n-1}]],

{v[kickcellN,kickcellM]'[t] ==(A*v[kickcellN,kickcellM][t]*(v[kickcellN,kickcellM][t]-α)*(1-v[kickcellN,kickcellM][t])-w[kickcellN,kickcellM][t]-w0)/ϵ+RC[[2 kickcellN -1, kickcellM-1]]*(K*(v[kickcellN,kickcellM-1][t]-v[kickcellN,kickcellM][t])) + RC[[2 kickcellN-2, kickcellM]]*(K*(v[kickcellN-1,kickcellM][t]-v[kickcellN,kickcellM][t]))+RC[[2 kickcellN-1,kickcellM]]*(K*(v[kickcellN,kickcellM+1][t]-v[kickcellN,kickcellM][t])) + RC[[2kickcellN,kickcellM]]*(K*(v[kickcellN+1,kickcellM][t]-v[kickcellN,kickcellM][t]))},

			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]-γpace*w[kickcellN,kickcellM][t]-v0pace},

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

(* Calculating pacemaker frequency *)

middlecell = Evaluate[v[kickcellN,kickcellM][t]/.fitz[[1]]];

middlecelltab = Table[Evaluate[middlecell],{t,0,T}];

cutoff = First[NMaximize[{v[kickcellN,kickcellM][t]/.fitz[[1]],0<t<T},t]]-0.2;

If[Count[Table[Abs[middlecelltab[[i]]-middlecelltab[[i-1]]]<10^-2,{i,190,T+1}],True]==0,

timeseries =Drop[ Union[Table[t/.FindRoot[middlecell==cutoff,{t,tInit,tmin,tmax}],{tInit,tmin+dt,tmax-dt,dt}],SameTest->(Abs[#1-#2]<10^-3&)],-1];

timeseries2 = Drop[timeseries, {2,Length[timeseries],2}];

timeperiod[JJ] = timeseries2[[Length[timeseries2]]]-timeseries2[[Length[timeseries2]-1]];,

timeperiod[JJ] = Null]];

periodtable2 = Table[timeperiod[JJ],{JJ,1,100,1}];

Export["Pacemaker_Frequency_K_"<>ToString[K]<>"_Prob_"<>ToString[prob]<>".xls",periodtable2]],JJ],prob]

(* Calculating a connectivity matrix *)

Conn = K*RC;

HConn =Transpose[Drop[Transpose[Table[Conn[[i]],{i,1,2n-1,2}]],-1]];
VConn = Table[Conn[[i]],{i,2,2n-1,2}];

HConn2 = HConn;
VConn2 = VConn;

For[i=n-1,i≥1,i--,
HConn2=Transpose[Insert[Transpose[HConn2],Transpose[HConn][[i]],i]]];

For[i=n-1,i≥1,i--,
VConn2=Insert[VConn2,VConn[[i]],i]];

(* Horizontal Connections *)

HConnCoord = Partition[Flatten[Table[{i,2j-1},{i,1,n},{j,1,n}]],2];
HConnMatrix = Insert[HConn2,-1,HConnCoord];
HConnCoord2 = Table[{i},{i,1,n+1}];
ZeroVector = Flatten[Table[{0},{i,1,3n-2}]];
HConnMatrix2 = Insert[HConnMatrix,ZeroVector,HConnCoord2];
HConnCoord3 = Table[{2i+1},{i,1,n-1}];
HConnMatrix3 = Insert[HConnMatrix2,ZeroVector,HConnCoord3];
HConnCoord4 = {{1},{3n}};
ZeroVector2 = Flatten[Table[{0},{i,1,3n}]];
HConnMatrix4 = Transpose[Insert[Transpose[HConnMatrix3],ZeroVector2,{{1},{3n-1}}]];

(* Vertical Connections *)

VConnCoord = Table[{2i-1},{i,1,n}];
NegCoordValue = Table[-1,{i,1,n}];
VConnMatrix = Insert[VConn2, NegCoordValue, VConnCoord];
VConnCoord2 = Table[{i},{i,1,n+1}];
VConnMatrix2 = Transpose[Insert[Transpose[VConnMatrix],ZeroVector,VConnCoord2]];
VConnCoord3 = Table[{2i+1},{i,1,n-1}];
VConnMatrix3 =Transpose[Insert[Transpose[VConnMatrix2],ZeroVector,VConnCoord3]];
VConnMatrix4 = Insert[VConnMatrix3,ZeroVector2,{{1},{3n-1}}];

ConnectionsMatrix =HConnMatrix4+VConnMatrix4;

NegVector = Table[-1,{i,1,3n}];
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{1}->NegVector];
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{3n}->NegVector];

For[i=2,i≤3n,i=i+1,
ConnectionsMatrix=ReplacePart[ConnectionsMatrix,{i,1}->-1];
ConnectionsMatrix=ReplacePart[ConnectionsMatrix,{i,3n}->-1]];

For[i=3,i≤3n,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j},{j,1,n}]->-1]];

For[i=4,i≤3n+1,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j},{j,1,n}]->-1]];

For[i=3,i≤3n,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j+1},{j,1,n}]->-1]];

For[i=4,i≤3n+1,i=i+3,
ConnectionsMatrix = ReplacePart[ConnectionsMatrix,Table[{i,3j+1},{j,1,n}]->-1]];

ConnectionsMatrix = ReplacePart[ConnectionsMatrix,{ (3n+1)/2,(3n+1)/2}->-3];

TruncatedMatrix = ConnectionsMatrix[[3*(kickcellN-6);;3*(kickcellN+6),3*(kickcellM-6);;3*(kickcellM+6)]];

TruncatedConnectionMatrixImage= MatrixPlot[TruncatedMatrix,ColorRules->{-3->White,-2->Darker[Black],-1->Lighter[Gray],0->Lighter[Gray]},ColorFunction->(Blend[{White,Darker[Gray]},#]&),AspectRatio->1]

Export["/Users/Rachel/Documents/MOAC/PhD/Modelling/FHN Model/Global Network Diagrams/Pacemaker_DecayingFrequency_Prob_0.6.jpg",TruncatedConnectionMatrixImage]