<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Simple model of epidemic dynamics
Step1:  Setup a python function that specifies the dynamics
Step2:  The function SIR above takes three arguments, $U$, $t$, and $p$ that represent the states of the system, the time and the parameters, respectively. 
Outbreak condition
The condition 
\begin{equation}
\frac{\alpha}{\beta}x(t)>1 , \quad y>0
\end{equation}
defines a threshold for a full epidemic outbreak. An equivalent condition is 
\begin{equation}
x>\frac{\beta}{\alpha }, \quad y>0
\end{equation}
Therefore, with the parameters $(\alpha,\beta)$=(0.5,0.1), there will be an outbreak if the initial condition for $x(t)>1/5$ with $y>0$. 
Notice that the initial value for $z$ can be interpreted as the initial proportion of immune individuals within the population. 
The dynamics related to the oubreak condition can be studied by defining a variable $B(t) = x(t) \alpha/\beta$, called by some authors "effective reproductive number". If $x(t)\approx 1$, the corresponding $B(t)$ is called "basic reproductive number", or $R_o$.
Let's define a python dictionary containing parameters and initial conditions to perform simulations.
Step3:  Integrate numerically and plot the results
<ASSISTANT_TASK:>
Python Code:
#Import the necessary modules and perform the necessary tests
import scipy as sc
import pylab as gr
sc.test("all",verbose=0)
%matplotlib inline
Explanation: Simple model of epidemic dynamics: SIR
Prof. Marco Arieli Herrera-Valdez,
Facultad de Ciencias, Universidad Nacional Autónoma de México
Created March 7, 2016
Let $x$, $y$, and $z$ represent the fraction of susceptibles, infected, and recovered individuals within a population. Assume homogeneous mixing with a probability of infection given a contact with an infected individual given by $\beta$ and an average removal time $\beta^{-1}$ from the infected group, by recovery or death due to infection. The population dynamics are given by
\begin{eqnarray}
\partial_t x &=& -\alpha xy
\
\partial_t y &=& \left( \alpha x - \beta \right) y
\
\partial_t x &=& \beta y
\end{eqnarray}
Notice that the population size does not matter because it is kept constant.
End of explanation
def SIR(U,t,p):
    x,y,z=U
    yNew= p["alpha"] * y * x
    zNew= p["beta"] * y 
    dx = -yNew
    dy = yNew - zNew
    dz = zNew
    return dx, dy, dz
Explanation: Setup a python function that specifies the dynamics
End of explanation
p={"alpha": 0.15, "beta":0.1, "timeStop":300.0, "timeStep":0.01 }
p["Ro"]=p["alpha"]/p["beta"]
p["sampTimes"]= sc.arange(0,p["timeStop"],p["timeStep"])
N= 1e4; i0= 1e1; r0=0; s0=N-i0-r0
x0=s0/N; y0=i0/N; z0=r0/N;
p["ic"]=[x0,y0,z0]
print("N=%g with initial conditions (S,I,R)=(%g,%g,%g)"%(N,s0,i0,r0))
print("Initial conditions: ", p["ic"])
print("B(0)=%g"%(p["ic"][0]*p["Ro"]))
Explanation: The function SIR above takes three arguments, $U$, $t$, and $p$ that represent the states of the system, the time and the parameters, respectively. 
Outbreak condition
The condition 
\begin{equation}
\frac{\alpha}{\beta}x(t)>1 , \quad y>0
\end{equation}
defines a threshold for a full epidemic outbreak. An equivalent condition is 
\begin{equation}
x>\frac{\beta}{\alpha }, \quad y>0
\end{equation}
Therefore, with the parameters $(\alpha,\beta)$=(0.5,0.1), there will be an outbreak if the initial condition for $x(t)>1/5$ with $y>0$. 
Notice that the initial value for $z$ can be interpreted as the initial proportion of immune individuals within the population. 
The dynamics related to the oubreak condition can be studied by defining a variable $B(t) = x(t) \alpha/\beta$, called by some authors "effective reproductive number". If $x(t)\approx 1$, the corresponding $B(t)$ is called "basic reproductive number", or $R_o$.
Let's define a python dictionary containing parameters and initial conditions to perform simulations.
End of explanation
# Numerical integration
xyz= sc.integrate.odeint(SIR, p["ic"], p["sampTimes"], args=(p,)).transpose()
# Calculate the outbreak indicator
B= xyz[0]*p["alpha"]/p["beta"]
# Figure
fig=gr.figure(figsize=(11,5))
gr.ioff()
rows=1; cols=2
ax=list()
for n in sc.arange(rows*cols):
    ax.append(fig.add_subplot(rows,cols,n+1))
ax[0].plot(p["sampTimes"], xyz[0], 'k', label=r"$(t,x(t))$")
ax[0].plot(p["sampTimes"], xyz[1], 'g', lw=3, label=r"$(t,y(t))$")
ax[0].plot(p["sampTimes"], xyz[2], 'b', label=r"$(t,z(t))$")
ax[0].plot(p["sampTimes"], B, 'r', label=r"$(t,B(t))$")
ax[0].plot([0, p["timeStop"]], [1,1], 'k--', alpha=0.4)
ax[1].plot(xyz[0], xyz[1], 'g', lw=3, label=r"$(x(t),y(t))$")
ax[1].plot(xyz[0], xyz[2], 'b', label=r"$(x(t),z(t))$")
ax[1].plot(xyz[0], B, 'r', label=r"$(x(t),B(t))$")
ax[1].plot([0, 1], [1,1], 'k--', alpha=0.4)
ax[0].legend(); ax[1].legend(loc="upper left")
gr.ion(); gr.draw()
Explanation: Integrate numerically and plot the results
End of explanation
<END_TASK>