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

(* Parameter values *)

K = 0.5;

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;

γpace = 0.1;
v0pace = 0.7;

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

(* Initiating variables *)

CC = Table[i,{i,1,1}];

inactivecount = Table[i,{i,1,1}];

prob = 0;
JJ = 1;
LL = 1;

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

(* Solving the FHN equations *)

Timing[

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

(* Counting the number of cells that are activated *)

CC[[JJ]] = MaxValue[{w[kickcellN,kickcellM][t]/.fitz[[1]],T-25≤t≤T},t];


activetable=CC-localmax;

(* Counting the number of inactive cells *)

inactivecount[[LL]]=Count[Negative[activetable],True]]