Below, you find the R code for computing the Game of life
n=20 # Size of matrix
mat=matrix(0,ncol=n,nrow=n) # Create a n x n matrix with zeros
mat[5:14,10]=1 # Add 10 live cells
temp_mat=mat # Create a temporary matrix
image(t(apply(mat, 2, rev)),col=c("grey50","seagreen1"),yaxt="n",xaxt="n") # Plot image
grid(nx=n,ny=n,col="grey70",lty=1)
mat1=mat[c(2:n,1),];mat2=mat[c(n,1:n-1),]; mat3=mat[,c(2:n,1)]; mat4=mat[,c(n,1:n-1)]; mat5=mat[c(2:n,1),c(2:n,1)]; mat6=mat[c(2:n,1),c(n,1:n-1)]
mat7=mat[c(n,1:n-1),c(n,1:n-1)];mat8=mat[c(n,1:n-1),c(2:n,1)];
for (k in 1:200){ # Repeat 200 times
numb_alive=mat1+mat2+mat3+mat4+mat5+mat6+mat7+mat8
temp_mat[mat==1 & numb_alive<2]=0
temp_mat[mat==1 & numb_alive>3]=0
temp_mat[mat==1 & (numb_alive==2 | numb_alive==3)]=1
temp_mat[mat==0 & numb_alive==3]=1
mat=temp_mat # Update matrix
mat1=mat[c(2:n,1),];mat2=mat[c(n,1:n-1),];mat3=mat[,c(2:n,1)];mat4=mat[,c(n,1:n-1)]
mat5=mat[c(2:n,1),c(2:n,1)];mat6=mat[c(2:n,1),c(n,1:n-1)];mat7=mat[c(n,1:n-1),c(2:n,1)]
mat8=mat[c(n,1:n-1),c(n,1:n-1)]
image(t(apply(mat, 2, rev)),col=c("grey50","seagreen1"),add=TRUE) # Plot image
grid(nx=n,ny=n,col="grey70",lty=1)
Sys.sleep(0.5) # To see changes on the screen we need to pause the loop
}
Below, you find the R code for the SIR model
rm(list=ls())
q=0.5 # Probability for recovery
n=50 # Size of matrix 50 x 50
mat=matrix(ncol=n,nrow=n,0) # Create matrix with zeros (Healthy individuals)
mat[25,15]=1 # Place an infected individual at pos (25,15)
image(mat,col=c("grey50","deeppink","seagreen1"),yaxt="n",xaxt="n",zlim=c(0,2)) # Plot image
grid(nx=n,ny=n,col="grey70",lty=1)
temp=mat
t=80 # Number of time steps
n_helalthy=rep(0,t) # Vector to count healthy at each time step
n_infected=rep(0,t) # Vector to count infected at each time step
n_resistant=rep(0,t) # Vector to count resistant at each time step
for (k in 1:t){ # Repeat t times
kH=0 # Initialize counter for healthy
kI=0 # Initialize counter for infected
kR=0 # Initialize counter for resistant
# Step through each element in the matrix
for (i in 1:n){
for (j in 1:n){
if (mat[i,j]==0) { kH=kH+1} # Count healthy
if (mat[i,j]==1) { kI=kI+1} # Count infected
if (mat[i,j]==2) { kR=kR+1} # Count resistant
R=0 # Initialize counter for number of infected neighbors
if (mat[i,j]==0){ # If healthy individual
E=i+1
W=i-1
N=j-1
S=j+1
# Check if outside the matrix
if (E==n+1) { E=1}
if (W==0) { W=n}
if (N==0) { N=n}
if (S==n+1) { S=1}
# Count number of infected neighbors
if(mat[E,j]==1){R=R+1} # East
if(mat[W,j]==1){R=R+1} # West
if(mat[i,N]==1){R=R+1} # North
if(mat[i,S]==1){R=R+1} # South
if(mat[E,N]==1){R=R+1} # North East
if(mat[E,S]==1){R=R+1} # South East
if(mat[W,N]==1){R=R+1} # North West
if(mat[W,S]==1){R=R+1} # South West
}
a=-1.5
b=0.6
Pinfect=(1/(1+exp(-(a+b*R)))) # Calc probability for healthy to become infected
g=runif(1) # Draw a random number between 0 and 1
if (g<Pinfect & mat[i,j]==0 & R>0){
temp[i,j]=1 # Healthy becomes infected
}
if (mat[i,j]==1){ # If infected individual
g=runif(1) # Draw a random number between 0 and 1
if (g<q){
temp[i,j]=2 # Infected becomes Resistant
}
}
}
}
image(mat,col=c("grey50","deeppink","seagreen1"),add=TRUE,,zlim=c(0,2)) # Plot image
grid(nx=n,ny=n,col="grey70",lty=1)
Sys.sleep(0.1) # To see movement on screen we need to pause the loop
mat=temp # Overwrite matrix
# Save number of healthy, infected and resistant at each time step
n_helalthy[k]=kH
n_infected[k]=kI
n_resistant[k]=kR
}
graphics.off()
plot(1:k,n_helalthy,type="l",ylab="Number",xlab="Time steps (weeks)",col=1,ylim=c(0,2600))
lines(1:k,n_infected,col=2)
lines(1:k,n_resistant,col=3)
legend(x=52,y=1599,c("Susceptible","Infected","Recovered"),lty=1,col=1:3)
Here is the code for the HIV model
rm(list=ls())
phiv=0.01 # Fraction of cells that are initially infected at time zero
Tao=5 # Number of time-steps a dead cell stays infected
p_rep = 0.97 # Probability that a dead cell is replaced by a healthy cell
n=100 # Size of matrix
NN=n*n # Total number of cells in grid
Time_A1=matrix(ncol=n,nrow=n,0) # A matrix that is used to count life time of A1 cells
mat=matrix(ncol=n,nrow=n,0) # Create matrix with zeros (Healthy CD4 T cells)
image(mat,col=c("grey50","yellow","red"),yaxt="n",xaxt="n",zlim=c(0,2)) # Plot image
# Add about phiv% (1%) of infected cells at random places
for (i in 1:n){
for (j in 1:n){
if (runif(1)<phiv) {mat[i,j]=1} # Add infected cell
}
}
image(mat,col=c("grey50","yellow","red"),yaxt="n",xaxt="n",zlim=c(0,2)) # Plot image
t=50 # Number of time-steps
temp=mat # Create temp matrix for synchronizing updating of matrix
n_healthy=rep(0,t) # Vector to count healthy cells at each time-step
n_infect=rep(0,t) # Vector to count infected cells at each time-step
n_dead=rep(0,t) # Vector to count dead cells at each time-step
n_healthy[1]=sum(mat==0) # Number of healthy cells at first time point
n_infect[1]=sum(mat==1) # Number of infected cells at first time point
n_dead[1]=sum(mat==2) # Number of dead cells at first time point
for (k in 1:t){ # Repeat t times
for (i in 1:n){
for (j in 1:n){
if (mat[i,j]==0){ # If healthy cell
E=i+1
W=i-1
N=j-1
S=j+1
# Check if outside the matrix
if (E==n+1) {E=1}
if (W==0) {W=n}
if (N==0) {N=n}
if (S==n+1) {S=1}
R=0 # Set to zero, number of surrounding infected cells
# Count number of infected neighbors
if(mat[E,j]==1){R=R+1} # East
if(mat[W,j]==1){R=R+1} # West
if(mat[i,N]==1){R=R+1} # North
if(mat[i,S]==1){R=R+1} # South
if(mat[E,N]==1){R=R+1} # North East
if(mat[E,S]==1){R=R+1} # South East
if(mat[W,N]==1){R=R+1} # North West
if(mat[W,S]==1){R=R+1} # South West
}
# Count how long infected has existed
if (mat[i,j]==1){ # If infected cell
Time_A1[i,j]=Time_A1[i,j]+1 # Count how long A1 has existed
}
# Rule 1: Update of infected cell
if (mat[i,j]==1 & Time_A1[i,j]>=Tao){ # If infected cell has existed for more than Tao time-steps
temp[i,j]=2 #infected -> Dead
Time_A1[i,j]=0
}
# Rule 2: Update of healthy cells
if (R>0 & mat[i,j]==0){
temp[i,j]=1 # Healthy becomes infected
}
### Rule 3: Dead cell is replaced by a Healthy cell or an infected cell ###
pr=runif(1)
if (mat[i,j]==2 & pr<p_rep){
temp[i,j]=0 # Dead -> Healthy
}
if (mat[i,j]==2 & pr>p_rep){
temp[i,j]=1 # Dead -> Infected
}
########################################################
}}
Sys.sleep(0.1) # To see movement on screen we need to pause the loop
image(mat,col=c("grey50","yellow","red"),yaxt="n",xaxt="n",zlim=c(0,2),add=TRUE) # Plot image
mat=temp # Overwrite matrix
# Save number of healthy, infected and dead cells at each time-step
n_healthy[k+1]=sum(mat==0) # count healthy cells
n_infect[k+1]=sum(mat==1) # count infected cells
n_dead[k+1]=sum(mat==2) # count dead cells
}
graphics.off()
plot(1:(t+1),n_healthy/NN,type="b",ylab="Cell Densities",xlab="time-steps",pch=0, lty=1,ylim=c(0,1),col="gray")
points(1:(t+1), n_infect/NN, type="b", pch=19, lty=1,col="yellow")
points(1:(t+1), n_dead/NN, type="b", pch=6, lty=1,col="red")
legend("topright",c("Healthy","Infected","Dead"),pch=c(0,19,6),col=c("gray","yellow","red"))