import numpy as np import sympy as sym import matplotlib.pyplot as plt from scipy.optimize import root from scipy.linalg import eig

#standard values of parameters N = 1.0 b = 0.018 alpha = 0.33 gamma = 0.2 e = 0.2 lamb = 0.1 beta0 = 1.68 mu = 0.1 kappa = 1117 #equations def beta(P): return beta0*(1-mu)*np.exp(-kappa*P/N) def Sdot(S, I, P): return b*(N-S) - (beta(P)*S*I)/N def Edot(S, E, I, P): return beta(P)*S*I/N - (alpha + b)*E def Idot(E, I): return alpha*E - (gamma + b)*I def Rdot(I, R): return gamma*I - b*R def Pdot(I, P): return e*gamma*I - lamb*P

# Runge Kutta 4 algorithm def RK4(S0, E0, I0, R0, P0, n, h=0.01): #n = amount of integration steps #h = size of time step assert abs(S0+E0+I0+R0 - N) < 1e-8, 'S0+E0+I0+R0 should be equal to N' #check if the starting conditions are valid #prepare arrays that will be returned and set first values equal to starting conditions S = np.zeros(n) S[0] = S0 E = np.zeros(n) E[0] = E0 I = np.zeros(n) I[0] = I0 R = np.zeros(n) R[0] = R0 P = np.zeros(n) P[0] = P0 #prepare arrays for derivatives in the calculated points and set first values equal to derivatives in the starting conditions dS = np.zeros(n) dS[0] = Sdot(S0, I0, P0) dE = np.zeros(n) dE[0] = Edot(S0, E0, I0, P0) dI = np.zeros(n) dI[0] = Idot(E0, I0) dR = np.zeros(n) dR[0] = Rdot(I0, R0) dP = np.zeros(n) dP[0] = Pdot(I0, P0) #prepare time array t = np.zeros(n) #runge kutta algorithm for i in range(1, n): k1 = h*Sdot(S[i-1], I[i-1], P[i-1]) l1 = h*Edot(S[i-1], E[i-1], I[i-1], P[i-1]) m1 = h*Idot(E[i-1], I[i-1]) n1 = h*Rdot(I[i-1], R[i-1]) o1 = h*Pdot(I[i-1], P[i-1]) k2 = h*Sdot(S[i-1]+k1/2, I[i-1]+m1/2, P[i-1]+o1/2) l2 = h*Edot(S[i-1]+k1/2, E[i-1]+l1/2, I[i-1]+m1/2, P[i-1]+o1/2) m2 = h*Idot(E[i-1]+l1/2, I[i-1]+m1/2) n2 = h*Rdot(I[i-1]+m1/2, R[i-1]+n1/2) o2 = h*Pdot(I[i-1]+m1/2, P[i-1]+o1/2) k3 = h*Sdot(S[i-1]+k2/2, I[i-1]+m2/2, P[i-1]+o2/2) l3 = h*Edot(S[i-1]+k2/2, E[i-1]+l2/2, I[i-1]+m2/2, P[i-1]+o2/2) m3 = h*Idot(E[i-1]+l2/2, I[i-1]+m2/2) n3 = h*Rdot(I[i-1]+m2/2, R[i-1]+n2/2) o3 = h*Pdot(I[i-1]+m2/2, P[i-1]+o2/2) k4 = h*Sdot(S[i-1]+k3, I[i-1]+m3, P[i-1]+o3) l4 = h*Edot(S[i-1]+k3, E[i-1]+l3, I[i-1]+m3, P[i-1]+o3) m4 = h*Idot(E[i-1]+l3, I[i-1]+m3) n4 = h*Rdot(I[i-1]+m3, R[i-1]+n3) o4 = h*Pdot(I[i-1]+m3, P[i-1]+o3) S[i] = S[i - 1] + (k1+(k2*2)+(k3*2)+k4)/6 E[i] = E[i - 1] + (l1+(l2*2)+(l3*2)+l4)/6 I[i] = I[i - 1] + (m1+(m2*2)+(m3*2)+m4)/6 R[i] = R[i - 1] + (n1+(n2*2)+(n3*2)+n4)/6 P[i] = P[i - 1] + (o1+(o2*2)+(o3*2)+o4)/6 t[i] = t[i - 1] + h dS[i] = Sdot(S[i], I[i], P[i]) dE[i] = Edot(S[i], E[i], I[i], P[i]) dI[i] = Idot(E[i], I[i]) dR[i] = Rdot(I[i], R[i]) dP[i] = Pdot(I[i], P[i]) return S, E, I, R, P, t, dS, dE, dI, dR, dP

#using the RK4 algorithm for certain starting conditions #evolution towards a stable point S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)

fig = plt.figure(figsize=(13,10)) ax1 = fig.add_subplot(321) ax1.plot(t, S, color='r') ax1.set_xlabel('t') ax1.set_ylabel('S') ax2 = fig.add_subplot(322) ax2.plot(t, E, color='g') ax2.set_xlabel('t') ax2.set_ylabel('E') ax3 = fig.add_subplot(323) ax3.plot(t, I, color='k') ax3.set_xlabel('t') ax3.set_ylabel('I') ax4 = fig.add_subplot(324) ax4.plot(t, R, color='y') ax4.set_xlabel('t') ax4.set_ylabel('R') ax5 = fig.add_subplot(325) ax5.plot(t, P, color='c') ax5.set_xlabel('t') ax5.set_ylabel('P') ax6 = fig.add_subplot(326) ax6.plot(t, S+E+I+R-N, color='m') ax6.set_xlabel('t') ax6.set_ylabel('S+E+I+R-N') plt.show()

#evolution towards a limit cycle alpha = 0.33 beta0 = 100 S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.82025, E0=0.0075, I0=0.01125, R0=0.1610, P0=0.0057, n=75000, h=0.01)

fig = plt.figure(figsize=(13,10)) ax1 = fig.add_subplot(321) ax1.plot(t, S, color='r') ax1.set_xlabel('t') ax1.set_ylabel('S') ax2 = fig.add_subplot(322) ax2.plot(t, E, color='g') ax2.set_xlabel('t') ax2.set_ylabel('E') ax3 = fig.add_subplot(323) ax3.plot(t, I, color='k') ax3.set_xlabel('t') ax3.set_ylabel('I') ax4 = fig.add_subplot(324) ax4.plot(t, R, color='y') ax4.set_xlabel('t') ax4.set_ylabel('R') ax5 = fig.add_subplot(325) ax5.plot(t, P, color='c') ax5.set_xlabel('t') ax5.set_ylabel('P') plt.show()

plt.plot(t, S+E+I+R-N) plt.title(r'Error on the condition $N=S+E+I+R$') plt.xlabel('t') plt.ylabel(r'S+E+I+R-N') plt.show()

alpha = 0.33 beta0 = 1.68 mu=0 S, E, I, R, P, t, dS, dE, dI, dR, dP = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01) mu=0.5 S1, E1, I1, R1, P1, t1, dS1, dE1, dI1, dR1, dP1 = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01) mu=1.0 S2, E2, I2, R2, P2, t2, dS2, dE2, dI2, dR2, dP2 = RK4(S0=0.9950, E0=0.0025, I0=0.0015, R0=0.0010, P0=1.0, n=75000, h=0.01)

plt.plot(t, I, label=r'$\mu=0$', color='r') plt.plot(t1, I1, label=r'$\mu=0.5$', color='g') plt.plot(t2, I2, label=r'$\mu=1$', color='b') plt.title('Influence of governmental measures') plt.xlabel('t') plt.ylabel('I') plt.legend() plt.show()

mu = 0.1 #reset mu

S, Ea, Ia, R, P, ta, dS, dE, dI, dR, dP = RK4(S0=0.99999, E0=0.00001, I0=0, R0=0, P0=0, n=25000, h=0.01) S, Eb, Ib, R, P, tb, dS, dE, dI, dR, dP = RK4(S0=0.95, E0=0, I0=0.05, R0=0.0, P0=0.0, n=25000, h=0.01) #S, Ec, Ic, R, P, tc, dS, dE, dI, dR, dP = RK4(S0=1.0, E0=-0.01, I0=0.01, R0=0.0, P0=0.0, n=75000, h=0.01) #S, Ed, Id, R, P, td, dS, dE, dI, dR, dP = RK4(S0=1.0, E0=0.01, I0=-0.01, R0=0.0, P0=0.0, n=75000, h=0.01)

plt.plot(Ea, Ia) #plt.plot(Eb, Ib) #plt.plot(Ec, Ic) #plt.plot(Ed, Id) plt.xlim(-0.02, 0.02) plt.ylim(-0.02, 0.02) #plt.plot(0.02, 0, 'ko') plt.plot(0.00001, 0, 'ko') plt.show()

alpha = 0.33 beta0 = 70 S0, E0, I0, R0, P0, t0, dS0, dE0, dI0, dR0, dP0 = RK4(S0=0.8483, E0=0.0031, I0=0.0058, R0=0.1428, P0=0.0047, n=25000, h=0.01) beta0 = 55 S1, E1, I1, R1, P1, t1, dS1, dE1, dI1, dR1, dP1 = RK4(S0=0.8561, E0=0.0050, I0=0.0075, R0=0.1314, P0=0.0044, n=25000, h=0.01) beta0 = 60 S2, E2, I2, R2, P2, t2, dS2, dE2, dI2, dR2, dP2 = RK4(S0=0.8479, E0=0.0032, I0=0.0089, R0=0.1400, P0=0.0054, n=25000, h=0.01) beta0 = 50 S3, E3, I3, R3, P3, t3, dS3, dE3, dI3, dR3, dP3 = RK4(S0=0.8567, E0=0.0055, I0=0.0093, R0=0.1285, P0=0.0046, n=25000, h=0.01) beta0 = 75 S4, E4, I4, R4, P4, t4, dS4, dE4, dI4, dR4, dP4 = RK4(S0=0.8258, E0=0.0048, I0=0.0164, R0=0.1530, P0=0.0067, n=25000, h=0.01) beta0 = 80 S5, E5, I5, R5, P5, t5, dS5, dE5, dI5, dR5, dP5 = RK4(S0=0.8123, E0=0.0228, I0=0.0221, R0=0.1428, P0=0.0047, n=25000, h=0.01) beta0 = 85 S6, E6, I6, R6, P6, t6, dS6, dE6, dI6, dR6, dP6 = RK4(S0=0.8064, E0=0.0236, I0=0.0240, R0=0.1460, P0=0.0049, n=25000, h=0.01) beta0 = 90 S7, E7, I7, R7, P7, t7, dS7, dE7, dI7, dR7, dP7 = RK4(S0=0.8314, E0=0.0016, I0=0.0074, R0=0.1596, P0=0.0061, n=25000, h=0.01) beta0 = 95 S8, E8, I8, R8, P8, t8, dS8, dE8, dI8, dR8, dP8 = RK4(S0=0.8110, E0=0.0048, I0=0.0188, R0=0.1654, P0=0.0075, n=25000, h=0.01) beta0 = 100 S9, E9, I9, R9, P9, t9, dS9, dE9, dI9, dR9, dP9 = RK4(S0=0.7932, E0=0.0227, I0=0.0288, R0=0.1553, P0=0.0057, n=25000, h=0.01) #print final positions and use these as starting conditions until we stay on a limit cycle ''' print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S0[-1], E0[-1], I0[-1], R0[-1], P0[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S1[-1], E1[-1], I1[-1], R1[-1], P1[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S2[-1], E2[-1], I2[-1], R2[-1], P2[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S3[-1], E3[-1], I3[-1], R3[-1], P3[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S4[-1], E4[-1], I4[-1], R4[-1], P4[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S5[-1], E5[-1], I5[-1], R5[-1], P5[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S6[-1], E6[-1], I6[-1], R6[-1], P6[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S7[-1], E7[-1], I7[-1], R7[-1], P7[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S8[-1], E8[-1], I8[-1], R8[-1], P8[-1])) print('S0={:.4f}, E0={:.4f}, I0={:.4f}, R0={:.4f}, P0={:.4f}'.format(S9[-1], E9[-1], I9[-1], R9[-1], P9[-1])) ''' #determine mins and maxes to use in bifurcation diagram beta0_values = [70, 55, 62, 50, 75, 80, 85, 90, 95, 100] maxes_beta = [np.max(I0), np.max(I1), np.max(I2), np.max(I3), np.max(I4), np.max(I5), np.max(I6), np.max(I7), np.max(I8), np.max(I9)] mins_beta = [np.min(I0), np.min(I1), np.min(I2), np.min(I3), np.min(I4), np.min(I5), np.min(I6), np.min(I7), np.min(I8), np.min(I9)]

#extra integrations to show spiraling towards the limit cycle/stable point beta0 = 100 S89, E89, I89, R89, P89, t89, dS89, dE89, dI89, dR89, dP89 = RK4(S0=0.82025, E0=0.0075, I0=0.01125, R0=0.1610, P0=0.0057, n=75000, h=0.01) S91, E91, I91, R91, P91, t91, dS91, dE91, dI91, dR91, dP89 = RK4(S0=0.7475, E0=0.0427, I0=0.0488, R0=0.1610, P0=0.0057, n=75000, h=0.01) beta0 = 20 Ss, Es, Is, Rs, Ps, ts, dSs, dEs, dIs, dRs, dPs = RK4(S0=0.7932, E0=0.0227, I0=0.0288, R0=0.1553, P0=0.0057, n=75000, h=0.01)

plt.plot(E89, I89, color='b') plt.plot(E91, I91, color='r') plt.plot(0.0075, 0.01125, 'ko') plt.plot(0.0427, 0.0488, 'ko') plt.plot(E9, I9, color='k') plt.ylabel('I') plt.xlabel('E') plt.title(r'$\beta_0=100$') plt.show() plt.plot(Es, Is, color='b') plt.plot(0.0227, 0.0288, 'ko') plt.ylabel('I') plt.xlabel('E') plt.title(r'$\beta_0=20$') plt.show() plt.plot(E3, I3, color='b', label=r'$\beta_0=50$') plt.plot(E1, I1, color='r', label=r'$\beta_0=55$') plt.plot(E2, I2, color='k', label=r'$\beta_0=60$') plt.plot(E0, I0, color='g', label=r'$\beta_0=70$') plt.plot(E4, I4, color='m', label=r'$\beta_0=75$') plt.ylabel('I') plt.xlabel('E') plt.legend() plt.show() plt.plot(E5, I5, color='g', label=r'$\beta_0=80$') plt.plot(E6, I6, color='r', label=r'$\beta_0=85$') plt.plot(E7, I7, color='k', label=r'$\beta_0=90$') plt.plot(E8, I8, color='b', label=r'$\beta_0=95$') plt.plot(E9, I9, color='m', label=r'$\beta_0=100$') plt.ylabel('I') plt.xlabel('E') plt.legend() plt.show()

plt.plot(P5, I5) plt.show()

beta0 = 50 alpha = 0.20 S10, E10, I10, R10, P10, t10, dS10, dE10, dI10, dR10, dP10 = RK4(S0=0.853318, E0=0.011345, I0=0.009508, R0=0.125829, P0=0.004264, n=50000, h=0.005) alpha = 0.25 S11, E11, I11, R11, P11, t11, dS11, dE11, dI11, dR11, dP11 = RK4(S0=0.851354, E0=0.005885, I0=0.010236, R0=0.132525, P0=0.005148, n=50000, h=0.005) alpha = 0.30 S12, E12, I12, R12, P12, t12, dS12, dE12, dI12, dR12, dP12 = RK4(S0=0.849711, E0=0.006169, I0=0.012282, R0=0.131838, P0=0.005212, n=50000, h=0.005) alpha = 0.22 S13, E13, I13, R13, P13, t13, dS13, dE13, dI13, dR13, dP13 = RK4(S0=0.855633, E0=0.009476, I0=0.008403, R0=0.126488, P0=0.004150, n=50000, h=0.005) alpha = 0.27 S14, E14, I14, R14, P14, t14, dS14, dE14, dI14, dR14, dP14 = RK4(S0=0.855109, E0=0.005119, I0=0.008625, R0=0.131147, P0=0.004823, n=50000, h=0.005) alpha = 0.34 S15, E15, I15, R15, P15, t15, dS15, dE15, dI15, dR15, dP15 = RK4(S0=0.853360, E0=0.007786, I0=0.011758, R0=0.127096, P0=0.004554, n=50000, h=0.005) alpha = 0.24 S16, E16, I16, R16, P16, t16, dS16, dE16, dI16, dR16, dP16 = RK4(S0=0.856571, E0=0.008556, I0=0.008048, R0=0.126825, P0=0.004107, n=50000, h=0.005) alpha = 0.195 S17, E17, I17, R17, P17, t17, dS17, dE17, dI17, dR17, dP17 = RK4(S0=0.852537, E0=0.011212, I0=0.009925, R0=0.126326, P0=0.004423, n=50000, h=0.005) alpha = 0.28 S18, E18, I18, R18, P18, t18, dS18, dE18, dI18, dR18, dP18 = RK4(S0=0.845182, E0=0.014548, I0=0.014180, R0=0.126090, P0=0.004276, n=50000, h=0.005) alpha = 0.32 S19, E19, I19, R19, P19, t19, dS19, dE19, dI19, dR19, dP19 = RK4(S0=0.848053, E0=0.008059, I0=0.014062, R0=0.129826, P0=0.005016, n=50000, h=0.005) #print final positions and use these as starting conditions until we stay on a limit cycle ''' print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S10[-1], E10[-1], I10[-1], R10[-1], P10[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S11[-1], E11[-1], I11[-1], R11[-1], P11[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S12[-1], E12[-1], I12[-1], R12[-1], P12[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S13[-1], E13[-1], I13[-1], R13[-1], P13[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S14[-1], E14[-1], I14[-1], R14[-1], P14[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S15[-1], E15[-1], I15[-1], R15[-1], P15[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S16[-1], E16[-1], I16[-1], R16[-1], P16[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S17[-1], E17[-1], I17[-1], R17[-1], P17[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S18[-1], E18[-1], I18[-1], R18[-1], P18[-1])) print('S0={:.6f}, E0={:.6f}, I0={:.6f}, R0={:.6f}, P0={:.6f}'.format(S19[-1], E19[-1], I19[-1], R19[-1], P19[-1])) ''' #determine mins and maxes to use in bifuration diagram alpha_values = [0.20, 0.25, 0.30, 0.22, 0.27, 0.34, 0.24, 0.19, 0.28, 0.32] maxes_alpha = [np.max(I10), np.max(I11), np.max(I12), np.max(I13), np.max(I14), np.max(I15), np.max(I16), np.max(I17), np.max(I18), np.max(I19)] mins_alpha = [np.min(I10), np.min(I11), np.min(I12), np.min(I13), np.min(I14), np.min(I15), np.min(I16), np.min(I17), np.min(I18), np.min(I19)]

beta0 = 50 alpha = 0.22 S89, E89, I89, R89, P89, t89, dS89, dE89, dI89, dR89, dP89 = RK4(S0=0.846612, E0=0.0115, I0=0.01125, R0=0.130638, P0=0.004150, n=300000, h=0.01) S91, E91, I91, R91, P91, t91, dS91, dE91, dI91, dR91, dP89 = RK4(S0=0.813562, E0=0.027, I0=0.0288, R0=0.130638, P0=0.004150, n=300000, h=0.01) alpha = 0.12 Ss, Es, Is, Rs, Ps, ts, dSs, dEs, dIs, dRs, dPs = RK4(S0=0.813562, E0=0.027, I0=0.0288, R0=0.130638, P0=0.004150, n=300000, h=0.01)

plt.plot(E89, I89, color='b') plt.plot(E91, I91, color='r') plt.plot(E13, I13, color='k') plt.plot(0.0115, 0.01125, 'ko') plt.plot(0.027, 0.0288, 'ko') plt.ylabel('I') plt.xlabel('E') plt.title(r'$\alpha=0.220$') plt.show() plt.plot(Es, Is, color='b') plt.plot(0.027, 0.0288, 'ko') plt.ylabel('I') plt.xlabel('E') plt.title(r'$\alpha=0.120$') plt.show() plt.plot(E17, I17, color='m', label=r'$\alpha=0.195$') plt.plot(E10, I10, color='b', label=r'$\alpha=0.200$') plt.plot(E13, I13, color='r', label=r'$\alpha=0.220$') plt.plot(E16, I16, color='k', label=r'$\alpha=0.240$') plt.plot(E11, I11, color='g', label=r'$\alpha=0.250$') plt.ylabel('I') plt.xlabel('E') plt.legend() plt.show() plt.plot(E14, I14, color='b', label=r'$\alpha=0.270$') plt.plot(E18, I18, color='r', label=r'$\alpha=0.280$') plt.plot(E12, I12, color='k', label=r'$\alpha=0.300$') plt.plot(E19, I19, color='g', label=r'$\alpha=0.320$') plt.plot(E15, I15, color='m', label=r'$\alpha=0.340$') plt.ylabel('I') plt.xlabel('E') plt.legend() plt.show()

def SEIRP_deriv(x): #x = [S,E,I,R,P] dS = b*(N-x[0]) - (beta(x[4])*x[0]*x[2])/N dE = beta(x[4])*x[0]*x[2]/N - (alpha + b)*x[1] dI = alpha*x[1] - (gamma + b)*x[2] dR = gamma*x[2] - b*x[3] dP = e*gamma*x[2] - lamb*x[4] return np.array([dS, dE, dI, dR, dP])

alpha=0.33 #fixed alpha method, options = 'hybr', {'maxfev': 100000} #choose method for finding roots start_point = (0.8567, 0.0055, 0.015, 0.1285, 0.0046) #define starting conditions betas = np.linspace(0, 100, 10000) #beta0 will be varied #prepare lists to save the solutions S_roots_beta = [] E_roots_beta = [] I_roots_beta = [] R_roots_beta = [] P_roots_beta = [] for beta0 in betas: #walk over all beta0-values sol = root(SEIRP_deriv, start_point, method=method, options=options) #find the roots if sol.success: #if a solution was found, add it to the lists S_roots_beta.append(sol.x[0]) E_roots_beta.append(sol.x[1]) I_roots_beta.append(sol.x[2]) R_roots_beta.append(sol.x[3]) P_roots_beta.append(sol.x[4]) else: print('solution not found')

plt.plot(betas[:4866], I_roots_beta[:4866], color='k') plt.plot(betas[4867:], I_roots_beta[4867:], 'k--') #The location of the Hopf-point, will be determined later on in the stability analysis plt.plot(beta0_values, maxes_beta, 'ko') plt.plot(beta0_values, mins_beta, 'ko') plt.xlim(0, 100) plt.ylim(-0.002, 0.025) plt.xlabel(r'$\beta_0$') plt.title(r'Bifurcation diagram for $\beta_0$') plt.ylabel('I') plt.plot(betas, np.zeros(betas.size), 'k--') #trivial solution (see paper) plt.show() plt.plot(betas, I_roots_beta, color='k') plt.xlim(0, 1) plt.ylim(-0.002, 0.004) plt.xlabel(r'$\beta_0$') plt.ylabel('I') plt.title(r'Zoom on the bifurcation diagram for $\beta_0$') plt.plot(betas, np.zeros(betas.size), 'k--') #trivial solution (see paper) plt.show()

beta0 = 50 #keep beta0 fixed method, options = 'hybr', {'maxfev': 100000} #choose method for finding roots start_point = (0.8000, 0.0050, 0.0050, 0.1850, 0.0050) #define starting conditions alphas = np.linspace(0, 0.5, 1000) #define range of alpha-values #prepare lists to save the solutions S_roots_alpha = [] E_roots_alpha = [] I_roots_alpha = [] R_roots_alpha = [] P_roots_alpha = [] for alpha in alphas: #walk over all alpha-values sol = root(SEIRP_deriv, start_point, method=method, options=options) #find the roots if sol.success: #if a solution was found, add it to the lists S_roots_alpha.append(sol.x[0]) E_roots_alpha.append(sol.x[1]) I_roots_alpha.append(sol.x[2]) R_roots_alpha.append(sol.x[3]) P_roots_alpha.append(sol.x[4]) else: print('solution not found')

plt.plot(alphas[:380], I_roots_alpha[:380], color='k') plt.plot(alphas[381:686], I_roots_alpha[381:686], 'k--') #The location of the first Hopf-point, will be determined later on in the stability analysis plt.plot(alphas[687:], I_roots_alpha[687:], color='k') #The location of the second Hopf-point, will be determined later on in the stability analysis plt.plot(alpha_values, maxes_alpha, 'ko') plt.plot(alpha_values, mins_alpha, 'ko') plt.xlim(0, 0.5) plt.ylim(-0.002, 0.017) plt.xlabel(r'$\alpha$') plt.title(r'Bifurcation diagram for $\alpha$') plt.ylabel('I') plt.plot(alphas, np.zeros(alphas.size), 'k--') #trivial solution (see paper) plt.show()

import sympy as sym #define SEIRP-model in sympy-format def SEIRP_sym(N, b, alpha, gamma, e, lamb, beta0, mu, kappa, S, E, I, R, P): eq1 = b*(N-S) - ((beta0*(1-mu)*sym.exp(-kappa*P/N))*S*I)/N eq2 = (beta0*(1-mu)*sym.exp(-kappa*P/N))*S*I/N - (alpha + b)*E eq3 = alpha*E - (gamma + b)*I eq4 = gamma*I - b*R eq5 = e*gamma*I - lamb*P return sym.Matrix([eq1, eq2, eq3, eq4, eq5]) #define symbols S_s = sym.Symbol('S') E_s = sym.Symbol('E') I_s = sym.Symbol('I') R_s = sym.Symbol('R') P_s = sym.Symbol('P') N_s = sym.Symbol('N') b_s = sym.Symbol('b') alpha_s = sym.Symbol('alpha') gamma_s = sym.Symbol('gamma') e_s = sym.Symbol('e') lamb_s = sym.Symbol('lambda') beta0_s = sym.Symbol('beta_0') mu_s = sym.Symbol('mu') kappa_s = sym.Symbol('kappa') Y = sym.Matrix([S_s, E_s, I_s, R_s, P_s]) #calculate Jacobian SEIRP_sym(N_s, b_s, alpha_s, gamma_s, e_s, lamb_s, beta0_s, mu_s, kappa_s, S_s, E_s, I_s, R_s, P_s).jacobian(Y)

#numpy-implementation of the Jacobian (easier to calculate eigenvalues) def jacobian(S, E, I, R, P): j = np.array([[-b-beta(P)*I/N, 0, -beta(P)*S/N, 0, I*S*beta(P)*kappa/N**2],[beta(P)*I/N, -alpha-b, beta(P)*S/N, 0, -I*S*beta(P)*kappa/N**2],[0, alpha, -gamma-b, 0, 0],[0, 0, gamma, -b, 0],[0, 0, e*gamma, 0, -lamb]]) return j

alpha = 0.33 #fixed alpha eigenvalues = np.zeros([np.size(betas), 5], dtype=float) #prepare an array to store the real parts of the eigenvalues for i, beta0 in enumerate(betas): #loop over beta0-values j = jacobian(N, 0, 0, 0, 0) #define jacobian for trivial solution eigs, _ = eig(j) #calculate eigenvalues with scipy (imaginary part is ignored because dtype of array is float) eigs = np.array(eigs) #make numpy array from tuple of eigenvalues eigenvalues[i] = np.real(eigs) #save in array

plt.plot(betas, eigenvalues[:, 0], 'b') plt.plot(betas, eigenvalues[:, 1], 'g') plt.plot(betas, eigenvalues[:, 2], 'r') plt.plot(betas, eigenvalues[:, 3], 'm') plt.plot(betas, eigenvalues[:, 4], 'k') plt.xlabel(r'$\beta_0$') plt.ylabel('eigenvalues') plt.title(r'Eigenvalues of the Jacobian when varrying $\beta_0$') plt.show() plt.plot(betas[0:100], eigenvalues[:, 0][0:100], 'b') plt.plot(betas[0:100], eigenvalues[:, 1][0:100], 'g') plt.plot(betas[0:100], eigenvalues[:, 2][0:100], 'r') plt.plot(betas[0:100], eigenvalues[:, 3][0:100], 'm') plt.plot(betas[0:100], eigenvalues[:, 4][0:100], 'k') plt.xlabel(r'$\beta_0$') plt.ylabel('eigenvalues') plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\beta_0$') plt.show() print('The black eigenvalues turns from negative to positive between beta0=0.25 and beta0=0.26 --> trivial solution is stable for beta0<=0.25') print('beta0={} --> black eigenvalue={}'.format(betas[25], eigenvalues[25, 4])) print('beta0={} --> black eigenvalue={}'.format(betas[26], eigenvalues[26, 4]))

beta0 = 50 #fixed beta0 eigenvalues = np.zeros([np.size(alphas), 5], dtype=float) #prepare an array to store the real parts of the eigenvalues for i, alpha in enumerate(alphas): #loop over alpha-values j = jacobian(N, 0, 0, 0, 0) #define jacobian for trivial solution eigs, _ = eig(j) #calculate eigenvalues with scipy (imaginary part is ignored because dtype of array is float) eigs = np.array(eigs) #make numpy array from tuple of eigenvalues eigenvalues[i] = np.real(eigs) #save in array

plt.plot(alphas, eigenvalues[:, 0], 'b') plt.plot(alphas, eigenvalues[:, 1], 'g') plt.plot(alphas, eigenvalues[:, 2], 'r') plt.plot(alphas, eigenvalues[:, 3], 'm') plt.plot(alphas, eigenvalues[:, 4], 'k') plt.xlabel(r'$\alpha$') plt.ylabel('eigenvalues') plt.title(r'Eigenvalues of the Jacobian when varrying $\alpha$') plt.show() plt.plot(alphas[0:100], eigenvalues[:, 0][0:100], 'b') plt.plot(alphas[0:100], eigenvalues[:, 1][0:100], 'g') plt.plot(alphas[0:100], eigenvalues[:, 2][0:100], 'r') plt.plot(alphas[0:100], eigenvalues[:, 3][0:100], 'm') plt.plot(alphas[0:100], eigenvalues[:, 4][0:100], 'k') plt.xlabel(r'$\alpha$') plt.ylabel('eigenvalues') plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\alpha$') plt.show() print('There is always a positive eigenvalue --> trivial solution is unstable for all alphas')

alpha = 0.33 #fixed alpha eigenvalues = np.zeros([np.size(betas), 5], dtype=float) #prepare array to store real parts of eigenvalues for i, beta0 in enumerate(betas): #loop over beta0-values and make indices available #use index to select a certain fixed point that was calulated earlier (see part on bifurcation diagrams) SS = S_roots_beta[i] EE = E_roots_beta[i] II = I_roots_beta[i] RR = R_roots_beta[i] PP = P_roots_beta[i] j = jacobian(SS, EE, II, RR, PP) #plug fixed point into jacobian eigs, _ = eig(j) #calculate eigenvalues (imaginary part is ignored because dtype of array is float) eigs = np.array(eigs) #make numpy array of tuple eigenvalues[i] = np.real(eigs) #store in array

plt.plot(betas, eigenvalues[:, 0], 'b') plt.plot(betas, eigenvalues[:, 1], 'g') plt.plot(betas, eigenvalues[:, 2], 'r') plt.plot(betas, eigenvalues[:, 3], 'm') plt.plot(betas, eigenvalues[:, 4], 'k') plt.xlabel(r'$\beta_0$') plt.ylabel('eigenvalues') plt.title(r'Eigenvalues of the Jacobian when varrying $\beta_0$') plt.show() plt.plot(betas, eigenvalues[:, 0], 'k', label='eigenvalue') plt.plot(betas, eigenvalues[:, 1], 'k') plt.plot(betas, eigenvalues[:, 2], 'k') plt.plot(betas, eigenvalues[:, 3], 'k') plt.plot(betas, eigenvalues[:, 4], 'k') plt.xlim(40, 100) plt.ylim(-0.01, 0.01) plt.hlines(0, 40, 100, label=r'$y=0$', color='r') plt.xlabel(r'$\beta_0$') plt.ylabel('eigenvalues') plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\beta_0$') plt.legend() plt.show()

print('eigenvalue = {}'.format(eigenvalues[:, 2][4866])) print('eigenvalue = {}'.format(eigenvalues[:, 2][4867])) print('beta_0 = {}'.format(betas[4866]))

beta0=50 #fixed beta0 eigenvalues = np.zeros([np.size(alphas), 5], dtype=float) #prepare array to store real parts of eigenvalues for i, alpha in enumerate(alphas): #loop over alpha-values and make indices available #use index to select a certain fixed point that was calulated earlier (see part on bifurcation diagrams) SS = S_roots_alpha[i] EE = E_roots_alpha[i] II = I_roots_alpha[i] RR = R_roots_alpha[i] PP = P_roots_alpha[i] j = jacobian(SS, EE, II, RR, PP) #plug fixed point into jacobian eigs, _ = eig(j) #calculate eigenvalues (imaginary part is ignored because dtype of array is float) eigs = np.array(eigs) #make numpy array of tuple eigenvalues[i] = np.real(eigs) #store in array

plt.plot(alphas, eigenvalues[:, 0], 'b') plt.plot(alphas, eigenvalues[:, 1], 'g') plt.plot(alphas, eigenvalues[:, 2], 'r') plt.plot(alphas, eigenvalues[:, 3], 'm') plt.plot(alphas, eigenvalues[:, 4], 'k') plt.xlabel(r'$\alpha$') plt.ylabel('eigenvalues') plt.title(r'Eigenvalues of the Jacobian when varrying $\alpha$') plt.show() plt.plot(alphas, eigenvalues[:, 0], 'k', label='eigenvalue') plt.plot(alphas, eigenvalues[:, 1], 'k') plt.plot(alphas, eigenvalues[:, 2], 'k') plt.plot(alphas, eigenvalues[:, 3], 'k') plt.plot(alphas, eigenvalues[:, 4], 'k') plt.xlim(0.18, 0.36) plt.ylim(-0.001, 0.001) plt.hlines(0, 0.18, 0.36, color='r', label=r'$y=0$') plt.legend() plt.xlabel(r'$\alpha$') plt.ylabel('eigenvalues') plt.title(r'Zoom on eigenvalues of the Jacobian when varrying $\alpha$') plt.show()

print('eigenvalue = {}'.format(eigenvalues[:, 2][380])) print('eigenvalue = {}'.format(eigenvalues[:, 2][381])) print('alpha = {}'.format(alphas[380])) print('\n') print('eigenvalue = {}'.format(eigenvalues[:, 2][686])) print('eigenvalue = {}'.format(eigenvalues[:, 2][687])) print('alpha = {}'.format(alphas[686]))

alpha = 0.33 beta0 = 70 S_array_grid = np.linspace(0, 0.0500, 25) E_array_grid = np.linspace(0, 0.0500, 25) I_array_grid = np.linspace(0, 0.0500, 25) R_array_grid = np.linspace(0, 0.0500, 25) P_array_grid = np.linspace(0, 0.0500, 25) B1 , C1 = np.meshgrid(E_array_grid, I_array_grid) #create grid DA1, DB1, DC1, DD1, DE1 = SEIRP_deriv([0.8304, B1, C1, 0.1497, 0.0065]) #change cross-section here M = (np.hypot(DB1, DC1)) # Norm of the growth rate M[ M == 0] = 0.0001 # Avoid zero division errors DB1 /= M # Normalize the arrows DC1 /= M #------------------------------------------------------- fig = plt.figure(figsize=(17,11)) plt.ylabel('I') plt.xlabel('E') Q = plt.quiver(B1, C1, DB1, DC1, M, pivot='middle', cmap=plt.cm.jet) plt.plot(E0,I0) plt.xlim(0, 0.0500) plt.ylim(0, 0.0500) plt.title('Velocity vector field in the E,I-plane, and a limit cycle that intersects this plane') plt.colorbar() plt.show()

plt.plot(E0, I0) for i in range(22000, 25000, 100): plt.quiver(E0[i], I0[i], dE0[i], dI0[i], pivot='tail', color='r') plt.xlabel('E') plt.ylabel('I') plt.title('Limit cycle with (E,I) velocity vectors') plt.show()