From 0604394eedebf5f4d601797d6e806af7d85fbae4 Mon Sep 17 00:00:00 2001 From: Andrew Guinness Date: Sun, 5 Nov 2017 22:06:39 -0500 Subject: [PATCH 1/2] Tut10 add --- Tut10.py | 102 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 Tut10.py diff --git a/Tut10.py b/Tut10.py new file mode 100644 index 0000000..7e6c5b2 --- /dev/null +++ b/Tut10.py @@ -0,0 +1,102 @@ +# Question 1 # + +# Loading Packages +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * + +# Simulating Equation +def ddSim(y, t0, r, K): + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] + +# Setting up parameters for each r condition +rlist = [-0.1,0.1,0.4,0.8,1] +# N0 and time range +N0 = 10 +K = 100 +times = range(0,50) +# creating a dataframe to store model outputs with a column for each r value +store_rs=pandas.DataFrame({"time":times, "r1":0, "r2":0, "r3":0, "r4":0, "r5":0}) + +# setting up the model +for i in range(0,len(rlist)): + pars=(rlist[i],K) + sim = spint.odeint(func=ddSim, y0=N0, t=times, args=pars) + store_rs.iloc[:,i]=sim[:,0] + +#Simulating plot +p = ggplot(store_rs, aes(x='time', y='r1')) + geom_line() #r1 +p = p + geom_line(aes(x='time', y='r2'), colour = "red") #r2 +p = p + geom_line(aes(x='time', y='r3'), colour = "blue") #r3 +p = p + geom_line(aes(x='time', y='r4'), colour = "green") #r4 +p = p + geom_line(aes(x='time', y='r5'), colour = "orange") #r5 +print p + +# Varying Ks +Klist = [10,50,100] +r = 0.2 +N0 =1 +store_Ks=pandas.DataFrame({"time":times, "K1":0, "K2":0, "K3":0}) + +for i in range(0,len(Klist)): + pars=(r,Klist[i]) + sim = spint.odeint(func=ddSim, y0=N0, t=times, args=pars) + store_Ks.iloc[:,i]=sim[:,0] + +p2 = ggplot(store_Ks, aes(x='time', y='K3')) + geom_line() +p2 = p2 + geom_line(aes(x='time',y='K2'), colour = "blue") +p2 = p2 + geom_line(aes(x='time', y='K1'), colour = "red") +print p2 + +#Varying Ns +Nlist=[1,50,100] +r=0.1 +K=50 +store_Ns=pandas.DataFrame({"time":times, "N0_1":0, "N0_2":0, "N0_3":0}) + +for i in range(0,len(Nlist)): + pars=(r,K) + sim = spint.odeint(func=ddSim, y0=Nlist[i], t=times, args=pars) + store_Ns.iloc[:,i]=sim[:,0] + +p3 = ggplot(store_Ns, aes(x='time',y='N0_3')) + geom_line() +p3 = p3 + geom_line(aes(x='time',y='N0_2'), colour = "green") +p3 = p3 + geom_line(aes(x='time',y='N0_1'), colour = 'yellow') +print p3 + +# Question 2 # + +# Multifaceted equation +def SIR(y,t0,beta,gamma): + S = y[0] + I = y[1] + R = y[2] + + dSdt = -beta*I*S + dIdt = beta*I*S - gamma*I + dRdt = gamma*I + return [dSdt, dIdt, dRdt] + +betalist = [0.0005, 0.005, 0.0001, 0.00005,0.0001,0.0002,0.0001] +gammalist = [0.05,0.5,0.1,0.1,0.05,0.05,0.06] +SIRvalues = [999,1,0] +times = range(0,500) + +SIRframe = pandas.DataFrame({"time":times, "S":0, "I":0, "R":0}) +#renaming columns because dumb +cols = ["S", "I", "R", "time"] +SIRframe = SIRframe[cols] + +for i in range(0,7): + pars = (betalist[i],gammalist[i]) + sim = spint.odeint(func=SIR, y0=SIRvalues, t=times, args=pars) + for j in range(0,3): + SIRframe.iloc[:,j]=sim[:,j] + +p4 = ggplot(SIRframe, aes(x='time', y='S')) + geom_line() #susceptible +p4 = p4 + geom_line(aes(x='time', y='I'), colour = "red") #infected +p4 = p4 + geom_line(aes(x='time', y='R'), colour = "blue") #removed +print p4 \ No newline at end of file From b87df68a992cfec974a0402feb07b6fdea14d882 Mon Sep 17 00:00:00 2001 From: Andrew Guinness Date: Thu, 9 Nov 2017 20:33:31 -0500 Subject: [PATCH 2/2] Updated tutorial & readme --- README.md | 16 +++++++++++++++- Tut10.py | 33 +++++++++++++++++++++------------ 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 897565a..85f3943 100644 --- a/README.md +++ b/README.md @@ -1 +1,15 @@ -# Intro_Biocom_ND_319_Tutorial10 \ No newline at end of file +# Intro_Biocom_ND_319_Tutorial10 + +QUESTION 2 COMMENTARY + +Unsurprisingly, as R0 increases (in our scenario, ranging from 0.5 to 10), the maximum daily prevalence and incidence rise along with it, and in the highest scenario, our disease affected nearly 100% of the population. + +The math for this mechanism is relatively simple: + +The disease is transmitted at some rate beta which forces people from susceptible into affected by the interaction of their quantities. (i.e. I increases with beta*I*S and S DECREASES in tandem). + +Additionally, people are removed from the infected pool at some rate gamma (either by recovery or death). Thus, I = beta*I*S - gamma*I. + + + + diff --git a/Tut10.py b/Tut10.py index 7e6c5b2..115bff9 100644 --- a/Tut10.py +++ b/Tut10.py @@ -27,6 +27,7 @@ def ddSim(y, t0, r, K): sim = spint.odeint(func=ddSim, y0=N0, t=times, args=pars) store_rs.iloc[:,i]=sim[:,0] +cols = ["black", "red", "blue", "green", "orange"] #Simulating plot p = ggplot(store_rs, aes(x='time', y='r1')) + geom_line() #r1 p = p + geom_line(aes(x='time', y='r2'), colour = "red") #r2 @@ -82,21 +83,29 @@ def SIR(y,t0,beta,gamma): betalist = [0.0005, 0.005, 0.0001, 0.00005,0.0001,0.0002,0.0001] gammalist = [0.05,0.5,0.1,0.1,0.05,0.05,0.06] +#inital values for S, I, R SIRvalues = [999,1,0] times = range(0,500) -SIRframe = pandas.DataFrame({"time":times, "S":0, "I":0, "R":0}) -#renaming columns because dumb -cols = ["S", "I", "R", "time"] -SIRframe = SIRframe[cols] +#creating a dataframe to store outputs for 7 scenarios +SIRframe = pandas.DataFrame({"beta":betalist, "gamma":gammalist, "MDI":0, "MDP":0, "percent_affected":0, "R0":0}) + +#calculating the R0 +SIRframe.R0 = (SIRframe.beta * 1000)/SIRframe.gamma + +inc = range(0,500) #initializing incidence parameter +prev = range(0,500) #initilaizing prevalence parameter for i in range(0,7): pars = (betalist[i],gammalist[i]) - sim = spint.odeint(func=SIR, y0=SIRvalues, t=times, args=pars) - for j in range(0,3): - SIRframe.iloc[:,j]=sim[:,j] - -p4 = ggplot(SIRframe, aes(x='time', y='S')) + geom_line() #susceptible -p4 = p4 + geom_line(aes(x='time', y='I'), colour = "red") #infected -p4 = p4 + geom_line(aes(x='time', y='R'), colour = "blue") #removed -print p4 \ No newline at end of file + sim = spint.odeint(func=SIR, y0=SIRvalues, t=times, args=pars) #simulate model for each beta,gamma pair + for j in range(500): + inc[j] = sim[j,1] - sim[j-1,1] #calculate incidence for all time points + SIRframe.iloc[i,0] = max(inc) #select maximum incidence + prev[j] = sim[j,1] / (sim[j,0] + sim[j,1] + sim[j,2]) #calculate prevalence for all time points + SIRframe.iloc[i,1] = max(prev) #select maximum prevalence + SIRframe.iloc[i,5] = (sim[499,1] + sim[499,2]) / (sim[499,0] + sim[499,1] + sim[499,2]) #calculate final percent affected + + +print SIRframe +