seircovid2S:

na = 28 # no of age groups
incub = 0
inf_per = 5 # total duration of human infectiousness
dur_cross_prot = 364
nu = DT/dur_cross_prot # used to model period of cross-protection 

init S[1..na,1] = Nb*(surv[i-1]-surv[i])/death[i]*exp(-(eq_FOI1+eq_FOI2+eq_FOI3+eq_FOI4)*mean_age[i])
init S[1..na,2..3] = 0

init I1[1..na,1..3] = 0
init I2[1..na,1..3] = 0

init R1[1..na,1] = Nb*(surv[i-1]-surv[i])/death[i]*(1-exp(-eq_FOI1*mean_age[i]))*exp(-(eq_FOI2+eq_FOI3+eq_FOI4)*mean_age[i])
init R2[1..na,1] = Nb*(surv[i-1]-surv[i])/death[i]*(1-exp(-eq_FOI2*mean_age[i]))*exp(-(eq_FOI1+eq_FOI3+eq_FOI4)*mean_age[i])

init I21[1..na,1..3] = 0
init I12[1..na,1..3] = 0

init R12[1..na,1] = Nb*(surv[i-1]-surv[i])/death[i]*(1-exp(-eq_FOI1*mean_age[i]))*(1-exp(-eq_FOI2*mean_age[i]))*exp(-(eq_FOI3+eq_FOI4)*mean_age[i])

Ntotal[1..na,1..3] = S[i,j]+I1[i,j]+I2[i,j]+R1[i,j]+R2[i,j]+I21[i,j]+I12[i,j]+R12[i,j]
NT = arraysum(Ntotal[*])

Y1[1..na,1..3] = (phi1[j]*inf_1[i,j]+phi21[j]*inf_21[i,j])
Y2[1..na,1..3] = (phi2[j]*inf_2[i,j]+phi12[j]*inf_12[i,j])

Y1T = arraysum(Y1[*])/NT
Y2T = arraysum(Y2[*])/NT

init infectious1 = 1e-7
init infectious2 = 1e-7
next infectious1 = delay(Y1T,incub,0)+infectious1-delay(Y1T,incub+inf_per,0)+extInfRand
next infectious2 = delay(Y2T,incub,0)+infectious2-delay(Y2T,incub+inf_per,0)+extInfRand

# Mwt_FOI1 = DT*Beta_hm_1*Kappa*infectious1
# Mwt_FOI2 = DT*Beta_hm_2*Kappa*infectious2
# Rel_R01 = 0.85
# Rel_R02 = 1
# Rel_R03 = 0.81
# Rel_R04 = 0.82
# Beta_hm_1 = R0agescale*Rel_R01*R0_req/(kappa*kappa*Mwt*inf_per*Beta_mh_mean/(1+delta*eip)/delta)
# Beta_hm_2 = R0agescale*Rel_R02*R0_req/(kappa*kappa*Mwt*inf_per*Beta_mh_mean/(1+delta*eip)/delta)
# Beta_hm_3 = R0agescale*Rel_R03*R0_req/(kappa*kappa*Mwt*inf_per*Beta_mh_mean/(1+delta*eip)/delta)
# Beta_hm_4 = R0agescale*Rel_R04*R0_req/(kappa*kappa*Mwt*inf_per*Beta_mh_mean/(1+delta*eip)/delta)

delta = 0.1  # adult mos death rate
eip = 8
R0_1 = kappa*kappa*Mwt*Beta_hm_1*inf_per*Beta_mh_mean/(1+delta*eip)/delta/R0agescale # eqn correct for exp distrib eip
R0_2 = kappa*kappa*Mwt*Beta_hm_2*inf_per*Beta_mh_mean/(1+delta*eip)/delta/R0agescale
eq_FOI1 = R0_1/lifespan # in time units of years
eq_FOI2 = R0_2/lifespan
EQT = STARTTIME
FOI1 = if(TIME<EQT) then DT*eq_FOI1/YL else DT*Beta_mh*Kappa*(Mwt_I1+Mwb_I1*Wb_relinf1)/NT
FOI2 = if(TIME<EQT) then DT*eq_FOI2/YL else DT*Beta_mh*Kappa*(Mwt_I2+Mwb_I2*Wb_relinf2)/NT

Acrit=3
FOIagescale=0.5
FOIas[1..na]=if(i>Acrit) then FOIagescale else 1

FOI1a[1..na]=FOIas[i]*FOI1
FOI2a[1..na]=FOIas[i]*FOI2

# update
O_S[1..na,1..3] = (rho1[j]*FOI1a[i]+rho2[j]*FOI2a[i]+agert[i]+deathrt[i]) * (S[i,j])
inf_1[1..na,1..3] = (rho1[j]*FOI1a[i]/(rho1[j]*FOI1a[i]+rho2[j]*FOI2a[i]+agert[i]+deathrt[i])) * (O_S[i,j])
inf_2[1..na,1..3] = (rho2[j]*FOI2a[i]/(rho2[j]*FOI2a[i]+agert[i]+deathrt[i]))*(O_S[i,j]-inf_1[i,j])

age_S[1..na,1..3]=(agert[i]/(agert[i]+deathrt[i]))*(O_S[i,j]-inf_1[i,j]-inf_2[i,j])

n_S[1..na,1] = ((if(i=1) then births else vacc_noncov_S[i-1,j]*age_S[i-1,j]) +S[i,j]-O_S[i,j])
n_S[1..na,2] = ((if((i=1)or(no_vacc_effect_after_inf=0)) then 0 else vacc_noncov_S[i-1,j]*age_S[i-1,j]+(1-vacc_noncov_S[i-1,1])*age_S[i-1,1]) +S[i,j]-O_S[i,j])
n_S[1..na,3] = ((if((i=1)or(no_vacc_effect_after_inf=1)) then 0 else vacc_noncov_S[i-1,j]*age_S[i-1,j]+(1-vacc_noncov_S[i-1,1])*age_S[i-1,1]) +S[i,j]-O_S[i,j])

O_I1[1..na,1..3] = (nu+agert[i]+deathrt[i])*(I1[i,j])
recov1[1..na,1..3] = (nu/(nu+agert[i]+deathrt[i]))*(O_I1[i,j])
age_I1[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_I1[i,j]-recov1[i,j])
age_I1[0,1..3] = 0
n_I1[1..na,1..3] = vacc_noncov[i-1,j]*age_I1[i-1,j]+nveai[j]*inf_1[i,j]+(if(j=3) then 0 else ((1-vacc_noncov[i-1,3-j])*age_I1[i-1,3-j])+nnveai[j]*inf_1[i,3-j])+I1[i,j]-O_I1[i,j]

O_I2[1..na,1..3] = (nu+agert[i]+deathrt[i])*(I2[i,j])
recov2[1..na,1..3] = (nu/(nu+agert[i]+deathrt[i]))*(O_I2[i,j])
age_I2[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_I2[i,j]-recov2[i,j])
age_I2[0,1..3] = 0
n_I2[1..na,1..3] = vacc_noncov[i-1,j]*age_I2[i-1,j]+nveai[j]*inf_2[i,j]+(if(j=3) then 0 else ((1-vacc_noncov[i-1,3-j])*age_I2[i-1,3-j])+nnveai[j]*inf_2[i,3-j])+I2[i,j]-O_I2[i,j]

O_R1[1..na,1..3] = (rho12[j]*FOI2a[i]+agert[i]+deathrt[i])*(R1[i,j])
inf_12[1..na,1..3] = (rho12[j]*FOI2a[i]/(rho12[j]*FOI2a[i]+agert[i]+deathrt[i]))*(O_R1[i,j])
age_R1[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_R1[i,j]-inf_12[i,j])
age_R1[0,1..3] = 0
n_R1[1..na,1..3] = vacc_noncov[i-1,j]*age_R1[i-1,j]+(if(j=3) then 0 else (1-vacc_noncov[i-1,3-j])*age_R1[i-1,3-j])+recov1[i,j]+R1[i,j]-O_R1[i,j]

O_R2[1..na,1..3] = (rho21[j]*FOI1a[i]+agert[i]+deathrt[i])*(R2[i,j])
inf_21[1..na,1..3] = (rho21[j]*FOI1a[i]/(rho21[j]*FOI1a[i]+agert[i]+deathrt[i]))*(O_R2[i,j])
age_R2[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_R2[i,j]-inf_21[i,j])
age_R2[0,1..3] = 0
n_R2[1..na,1..3] = vacc_noncov[i-1,j]*age_R2[i-1,j]+(if(j=3) then 0 else (1-vacc_noncov[i-1,3-j])*age_R2[i-1,3-j])+recov2[i,j]+R2[i,j]-O_R2[i,j]

O_I12[1..na,1..3] = (nu+agert[i]+deathrt[i])*(I12[i,j])
recov12[1..na,1..3] = (nu/(nu+agert[i]+deathrt[i]))*(O_I12[i,j])
age_I12[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_I12[i,j]-recov12[i,j])
age_I12[0,1..3] = 0
n_I12[1..na,1..3] = vacc_noncov[i-1,j]*age_I12[i-1,j]+(if(j=3) then 0 else (1-vacc_noncov[i-1,3-j])*age_I12[i-1,3-j])+inf_12[i,j]+I12[i,j]-O_I12[i,j]

O_I21[1..na,1..3] = (nu+agert[i]+deathrt[i])*(I21[i,j])
recov21[1..na,1..3] = (nu/(nu+agert[i]+deathrt[i]))*(O_I21[i,j])
age_I21[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_I21[i,j]-recov21[i,j])
age_I21[0,1..3] = 0
n_I21[1..na,1..3] = vacc_noncov[i-1,j]*age_I21[i-1,j]+(if(j=3) then 0 else (1-vacc_noncov[i-1,3-j])*age_I21[i-1,3-j])+inf_21[i,j]+I21[i,j]-O_I21[i,j]

O_R12[1..na,1..3] = (agert[i]+deathrt[i])*(R12[i,j])
age_R12[1..na,1..3] = (agert[i]/(agert[i]+deathrt[i]))*(O_R12[i,j])
age_R12[0,1..3] = 0
n_R12[1..na,1..3] = vacc_noncov[i-1,j]*age_R12[i-1,j]+(if(j=3) then 0 else (1-vacc_noncov[i-1,3-j])*age_R12[i-1,3-j])+recov12[i,j]+recov21[i,j]+R12[i,j]-O_R12[i,j]

next S[1..na,1] = vcu_noncov[i,j]*n_S[i,j]
next S[1..na,2] = (if(no_vacc_effect_after_inf=0) then 0 else (1-vcu_noncov[i,1]))*n_S[i,1]+n_S[i,j]
next S[1..na,3] = (if(no_vacc_effect_after_inf=1) then 0 else (1-vcu_noncov[i,1]))*n_S[i,1]+n_S[i,j]

next I1[1..na,1] = vcu_noncov[i,j]*n_I1[i,j]
next I1[1..na,2] = (1-vcu_noncov[i,1])*n_I1[i,1]+n_I1[i,2]
next I1[1..na,3] = n_I1[i,j]

next I2[1..na,1] = vcu_noncov[i,j]*n_I2[i,j]
next I2[1..na,2] = (1-vcu_noncov[i,1])*n_I2[i,1]+n_I2[i,2]
next I2[1..na,3] = n_I2[i,j]

next R1[1..na,1] = vcu_noncov[i,j]*n_R1[i,j]
next R1[1..na,2] = (1-vcu_noncov[i,1])*n_R1[i,1]+n_R1[i,2]
next R1[1..na,3] = n_R1[i,j]

next R2[1..na,1] = vcu_noncov[i,j]*n_R2[i,j]
next R2[1..na,2] = (1-vcu_noncov[i,1])*n_R2[i,1]+n_R2[i,2]
next R2[1..na,3] = n_R2[i,j]

next I12[1..na,1] = vcu_noncov[i,j]*n_I12[i,j]
next I12[1..na,2] = (1-vcu_noncov[i,1])*n_I12[i,1]+n_I12[i,2]
next I12[1..na,3] = n_I12[i,j]

next I21[1..na,1] = vcu_noncov[i,j]*n_I21[i,j]
next I21[1..na,2] = (1-vcu_noncov[i,1])*n_I21[i,1]+n_I21[i,2]
next I21[1..na,3] = n_I21[i,j]

next R12[1..na,1] = vcu_noncov[i,j]*n_R12[i,j]
next R12[1..na,2] = (1-vcu_noncov[i,1])*n_R12[i,1]+n_R12[i,2]
next R12[1..na,3] = n_R12[i,j]