Toy SIR model to illustrate role of transmissibility and generation time on growth of competing variants of an infectious pathogen

growthExample.R

library( deSolve )
library( ggplot2 )
library( cowplot )

dy <- function( t, y, parms, ... )
{
	S = y[1] 
	I0 = y[2] 
	I1 = y[3] 
	R = y[4] 
	psi = y[5]
	with ( as.list( parms ),{
		list(c(
		S = - (beta0*I0+beta1*I1)*S/N
		, I0 = beta0*I0*S/N - gamma0*I0 
		, I1 = beta1*I1*S/N - gamma1*I1
		, R = gamma0*I0 + gamma1*I1
		, psi = (S/N)*(beta1-beta0) + (gamma0-gamma1)
		))
	})
}

s=.5
parms_a = list(
	N = 1
	, beta0=1.25
	, beta1 = (1+s)*1.25
	, gamma0=1
	, gamma1 = 1
)
r=.5 * 1.25/(1.25-1)
parms_b = list(
	N = 1
	, beta0=1.25
	, beta1 = (1+r)*1.25
	, gamma0=1
	, gamma1 = (1+r)*1
)

times <- seq(0, 18, length=50)
I0 = 1e-3
I1 = 1e-4
y0 <- c( S = 1-(I0+I1)
  , I0=I0
  , I1=I1
  , R = 0
  , psi= log(I1/I0) )
a = ode( y0, times, dy, parms_a )
#~ plot(a)
b = ode( y0, times, dy, parms_b )
#~ plot(b)



A = as.data.frame(a) |> data.frame( Scenario='Higher transmissibility')
A$Infections <- A$I0 + A$I1
B = as.data.frame(b) |> data.frame( Scenario='Shorter generation time')
B$Infections <- B$I0 + B$I1
X = rbind( A, B )

cols <- c(`Higher transmissibility` = '#f69543ff'
, `Shorter generation time` = '#009ea2ff' ) 
p1 <- ggplot( data=X, aes(x=time,y=Infections,colour=Scenario) ) + geom_path(linewidth=2) + theme_classic() + scale_y_log10() + xlab('' ) + scale_colour_manual( values = cols )
p2 <- ggplot( data=X, aes(x=time,y=psi,colour=Scenario) ) + geom_path(linewidth=2) + theme_classic() + xlab('Time')+ scale_colour_manual( values = cols ) + ylab(expression(paste(sep='',psi,'(t)'))) #ylab('\u1D6F9')
print(p1)
print(p2)

p1 <- p1  + theme(legend.position = 'none' )

p2 <- p2 + theme( legend.position = c(.3,.99) 
, legend.key = element_rect(fill = alpha("white", 0.0))
, legend.background = element_rect(fill = alpha("white", 0.0))
)

p12 = plot_grid( plotlist = list( p1, p2 ) , nrow=2
, labels = c( 'a.', 'b.' )
)
ggsave( p12, file = 'R.svg', width=3.5,height=5.5 ) 
ggsave( p12, file = 'R.pdf', width=3.5,height=5.5 ) 

print( p12 )