Examan 1
Alumno: César Alfaro
Pneumonic Plague
Instrucciones:
LEE CUIDADOSAMENTE:
Para cada una de las preguntas abiertas, responde sólo lo que se solicita.
El examen consta solamente de una parte práctica, con un valor de 100 puntos.
Lee cuidadosamente las instrucciones del caso de estudio. Recuerda que debes tener a la mano la plantilla que utilizamos para generar modelos en R, así como Vensim listo para trabajar.
Recuerda que deberás adjuntar imagenes del diagrama realizado en Vensim. Si no se integra todo en un solo archivo html la calificación final del examen se penalizará calificándose sobre el 50%.
Puedes consultar tus notas, laboratorios, tareas y el libro de texto.
¡Éxito!
Introducción:
On 3 August 2009, media reported about an outbreak of pneumonic plague in north-west China:
“A second man has died of pneumonic plague in a remote part of north-west China where a town of more than 10,000 people has been sealed off. […] Local officials in north-western China have told the BBC that the situation is under control, and that schools and offices are open as usual. But to prevent the plague [from] spreading, the authorities have sealed off Ziketan, which has some 10,000 residents. About 10 other people inside the town have so far contracted the disease, according to state media. No-one is being allowed [to] leave the area, and the authorities are trying to track down people who had contact with the men who died. […] According to the WHO, pneumonic plague is the most virulent and least common form of plague. It is caused by the same bacteria that occur in bubonic plague – the Black Death that killed an estimated 25 million people in Europe during the Middle Ages. But while bubonic plague is usually transmitted by flea bites and can be treated with antibiotics, [pneumonic plague, which attacks the lungs, can spread from person to person or from animals to people], is easier to contract and if untreated, has a very high case-fatality ratio.” (BBC 03/08/2009)
Descripción del caso:
You are asked to make a SD model of this outbreak. Use the following assumptions: The total population of Ziketan amounted initially to 10000 citizens. New infections make that citizens belonging to the susceptible population become part of the infected population, which initially consists of just 1 person. The number of infections equals the product of the infection ratio, the contact rate, the susceptible population, and the infected fraction. Initially, the normal contact rate amounts to 50 contacts per week and the infection ratio to a staggering 75% per contact. The infected fraction equals of course the infected population over the sum of all other subpopulations. If citizens from the infected population die, they enter the statistics of the deceased population, else they are quarantined to recover. The recovering could be modeled simplistically as (1 − fatality ratio) ∗ infected population / recovery time. Suppose for the sake of simplicity that the average recovery time and the average decease time are both 2 days. The fatality ratio depends on the antibiotics coverage of the population which –in this poor part of China– is 0% at first.
The fatality ratio is 90% at 0% antibiotics coverage of the population and 15% at 100% antibiotics coverage of the population. Assume for the sake of simplicity that those belonging to the recovering population do not pose any threat of infection, either because they are really quarantined or because they are not contagious any more.
Preguntas del caso:
Desarrolla el diagrama de flujo del caso (25 puntos, producto: diagrama de flujo)
Construye el modelo de dinámica de sistemas basado en R y simula el modelo por un periodo de 1 mes.(25 puntos, productos: hipótesis dinámica y gráfico que la respalde)
Grafica y comenta la evolución de infections, deaths, recovering population y deceased population (25 puntos, producto: gráficos de evolución y comentarios del comportamiento)
Suponga ahora que la antibiotics coverage of the population es del 100%. Compara la evolución de de infections, deaths, recovering population y deceased population con las obtenidas en la pregunta anterior. (25 puntos, productos: código en R, gráficas y comentarios comparando ambos modelos)
Solución
Diagrama de flujo
Modelo dinámico del sistema para Pneumonic Plague
library(deSolve)
# Condiciones iniciales
InitialConditions <- c(susceptible.population = 9999, # people
infected.population = 1, # people
deceased.population = 0, # people
recovering.population = 0) # people
# Vector de tiempo (1 mes = 30 días)
times <- seq(0, # tiempo inicial, días
30, # tiempo final, días
0.25) # paso de tiempo, días
# Función del modelo
pneumonic.plague <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Endogenous auxiliary variables
# The infected fraction represents the probability of contact with infected person
infected.fraction <- infected.population / total.population # dimensionless
# Fatality ratio depends on antibiotics coverage (linear interpolation)
# 90% at 0% coverage, 15% at 100% coverage
fatality.ratio <- 0.9 - 0.75 * antibiotics.coverage.population # dimensionless
# Flow variables
# New infections equal the product of infection ratio, contact rate,
# susceptible population, and infected fraction
infections <- infection.ratio * normal.contact.rate * susceptible.population * infected.fraction # people/day
# Recovering population flow (those who survive the infection)
recovering <- (1 - fatality.ratio) * infected.population / recovery.time # people/day
# Deaths flow (those who die from the infection)
deaths <- fatality.ratio * infected.population / decease.time # people/day
# State (stock) variables
# Susceptible population decreases due to new infections
dsusceptible.population <- (-1) * infections # people/day
# Infected population increases with new infections and decreases with recoveries and deaths
dinfected.population <- infections - recovering - deaths # people/day
# Deceased population increases with deaths
ddeceased.population <- deaths # people/day
# Recovering population increases with recoveries
drecovering.population <- recovering # people/day
# List of state variables for integration and additional outputs
list(c(dsusceptible.population,
dinfected.population,
ddeceased.population,
drecovering.population),
infections = infections,
deaths = deaths,
recovering = recovering)
})
}
# Parámetros del modelo
parameters <- c(total.population = 10000, # people
infection.ratio = 0.75, # dimensionless (75% per contact)
normal.contact.rate = 50/7, # contacts/day (50 contacts/week)
recovery.time = 2, # days
decease.time = 2, # days
antibiotics.coverage.population = 0) # dimensionless (0% coverage)
# Método de integración
intg.method <- c("rk4")
# Ejecutar simulación
out <- ode(y = InitialConditions,
times = times,
func = pneumonic.plague,
parms = parameters,
method = intg.method)
# Gráfico de resultados
plot(out, col = c("blue"))
# Hipótesis dinámica:
# El modelo representa un brote epidémico de peste neumónica en una población cerrada de 10,000 habitantes.
# Con una infectividad del 75% por contacto, una tasa de contacto de aproximadamente 7 contactos/día, y sin cobertura de antibióticos (tasa de mortalidad del 90%), se espera observar un crecimiento exponencial inicial de infecciones seguido de una alta mortalidad que limita la propagación. La población susceptible disminuirá rápidamente, mientras que la población fallecida crecerá significativamente debido a la alta mortalidaad.
- Evolución de las variables infections, deaths, recovering population y _deceased population
# Usar el mismo modelo de la pregunta 2
library(deSolve)
InitialConditions <- c(susceptible.population = 9999, # people
infected.population = 1, # people
deceased.population = 0, # people
recovering.population = 0) # people
times <- seq(0, 30, 0.25) # 30 días
pneumonic.plague <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Endogenous auxiliary variables
infected.fraction <- infected.population / total.population # dimensionless
fatality.ratio <- 0.9 - 0.75 * antibiotics.coverage.population # dimensionless
# Flow variables
infections <- infection.ratio * normal.contact.rate * susceptible.population * infected.fraction # people/day
recovering <- (1 - fatality.ratio) * infected.population / recovery.time # people/day
deaths <- fatality.ratio * infected.population / decease.time # people/day
# State (stock) variables
dsusceptible.population <- (-1) * infections # people/day
dinfected.population <- infections - recovering - deaths # people/day
ddeceased.population <- deaths # people/day
drecovering.population <- recovering # people/day
list(c(dsusceptible.population,
dinfected.population,
ddeceased.population,
drecovering.population),
infections = infections,
deaths = deaths,
recovering = recovering)
})
}
parameters <- c(total.population = 10000, # people
infection.ratio = 0.75, # dimensionless
normal.contact.rate = 50/7, # contacts/day
recovery.time = 2, # days
decease.time = 2, # days
antibiotics.coverage.population = 0) # dimensionless
intg.method <- c("rk4")
out <- ode(y = InitialConditions,
times = times,
func = pneumonic.plague,
parms = parameters,
method = intg.method)
# Gráficos específicos de las variables solicitadas
# Para escribir el código de las gráficas obtuve ayuda de Claude AI
# Gráfico 1: Variables de flujo (infections y deaths)
par(mfrow = c(2, 2)) # 2x2 grid de gráficos
plot(out[, "time"], out[, "infections"],
type = "l", col = "red", lwd = 2,
xlab = "Time (days)", ylab = "New Infections (people/day)",
main = "Evolution of New Infections")
grid()
plot(out[, "time"], out[, "deaths"],
type = "l", col = "black", lwd = 2,
xlab = "Time (days)", ylab = "Deaths (people/day)",
main = "Evolution of Deaths")
grid()
# Gráfico 2: Variables de estado (recovering population y deceased population)
plot(out[, "time"], out[, "recovering.population"],
type = "l", col = "green", lwd = 2,
xlab = "Time (days)", ylab = "Recovering Population (people)",
main = "Evolution of Recovering Population")
grid()
plot(out[, "time"], out[, "deceased.population"],
type = "l", col = "purple", lwd = 2,
xlab = "Time (days)", ylab = "Deceased Population (people)",
main = "Evolution of Deceased Population")
grid()
# Resetear layout
par(mfrow = c(1, 1))
# Comentarios sobre el comportamiento:
# INFECTIONS (Nuevas infecciones por día):
# La curva muestra un brote muy muy alto con un pico de alrededor de 11 mil personas/día en el día 2. Parece que esta tasa excede incluso a la población total debido a la alta infectividad (75%) y frecuencia de contacto. El pico indica que la epidemia se propaga instantáneamente al inicio, seguido de una caída abrupta conforme se reduce la población susceptible y la alta mortalidad elimina a los infectados.
# DEATHS (Muertes por día):
# El pico de muertes ocurre en el día 2.5 con alrededor de 3 mil personas/día, se observa algo retrasado respecto a las infecciones debido al tiempo de progresión de la enfermedad (2 días).
# La alta tasa de mortalidad (90%) sin antibióticos hace que mucha gente infectada se convierte en fallecimientos.
# RECOVERING POPULATION (Población recuperada acumulada):
# La población recuperada crece hasta alcanzar alrededor de mil personas al final del mes, representando exactamente el 10% de la población total.
# DECEASED POPULATION (Población fallecida acumulada):
# La curva sugiere que que se alcanzan alrededor de 9 mil personas fallecidas (90% de la
# población total).- Comparación del modelo anterior contra la opción den donde se supone que la varaible antibiotics coverage of the population es del 100%. Compara la evolución de de infections, deaths, recovering population y deceased population
library(deSolve)
InitialConditions <- c(susceptible.population = 9999, # people
infected.population = 1, # people
deceased.population = 0, # people
recovering.population = 0) # people
times <- seq(0, 30, 0.25) # 30 días
pneumonic.plague <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Endogenous auxiliary variables
infected.fraction <- infected.population / total.population # dimensionless
fatality.ratio <- 0.9 - 0.75 * antibiotics.coverage.population # dimensionless
# Flow variables
infections <- infection.ratio * normal.contact.rate * susceptible.population * infected.fraction # people/day
recovering <- (1 - fatality.ratio) * infected.population / recovery.time # people/day
deaths <- fatality.ratio * infected.population / decease.time # people/day
# State (stock) variables
dsusceptible.population <- (-1) * infections # people/day
dinfected.population <- infections - recovering - deaths # people/day
ddeceased.population <- deaths # people/day
drecovering.population <- recovering # people/day
list(c(dsusceptible.population,
dinfected.population,
ddeceased.population,
drecovering.population),
infections = infections,
deaths = deaths,
recovering = recovering)
})
}
# Parámetros para escenario SIN antibióticos (0% cobertura)
parameters_no_antibiotics <- c(total.population = 10000, # people
infection.ratio = 0.75, # dimensionless
normal.contact.rate = 50/7, # contacts/day
recovery.time = 2, # days
decease.time = 2, # days
antibiotics.coverage.population = 0) # 0% coverage
# Parámetros para escenario CON antibióticos (100% cobertura)
parameters_with_antibiotics <- c(total.population = 10000, # people
infection.ratio = 0.75, # dimensionless
normal.contact.rate = 50/7, # contacts/day
recovery.time = 2, # days
decease.time = 2, # days
antibiotics.coverage.population = 1) # 100% coverage
intg.method <- c("rk4")
# Simulación SIN antibióticos
out_no_antibiotics <- ode(y = InitialConditions,
times = times,
func = pneumonic.plague,
parms = parameters_no_antibiotics,
method = intg.method)
# Simulación CON antibióticos
out_with_antibiotics <- ode(y = InitialConditions,
times = times,
func = pneumonic.plague,
parms = parameters_with_antibiotics,
method = intg.method)
# Gráficas comparativas
par(mfrow = c(2, 2))
# En todos los casos para construir las gráficas comparativas use Claude ID para mejorar el código.
# 1. Comparación de INFECTIONS
plot(out_no_antibiotics[, "time"], out_no_antibiotics[, "infections"],
type = "l", col = "red", lwd = 2,
xlab = "Time (days)", ylab = "New Infections (people/day)",
main = "Infections: 0% vs 100% Antibiotics Coverage",
ylim = c(0, max(out_no_antibiotics[, "infections"])))
lines(out_with_antibiotics[, "time"], out_with_antibiotics[, "infections"],
col = "blue", lwd = 2)
legend("topright", legend = c("0% Antibiotics", "100% Antibiotics"),
col = c("red", "blue"), lwd = 2)
grid()
# 2. Comparación de DEATHS
plot(out_no_antibiotics[, "time"], out_no_antibiotics[, "deaths"],
type = "l", col = "red", lwd = 2,
xlab = "Time (days)", ylab = "Deaths (people/day)",
main = "Deaths: 0% vs 100% Antibiotics Coverage",
ylim = c(0, max(out_no_antibiotics[, "deaths"])))
lines(out_with_antibiotics[, "time"], out_with_antibiotics[, "deaths"],
col = "blue", lwd = 2)
legend("topright", legend = c("0% Antibiotics", "100% Antibiotics"),
col = c("red", "blue"), lwd = 2)
grid()
# 3. Comparación de RECOVERING POPULATION
plot(out_no_antibiotics[, "time"], out_no_antibiotics[, "recovering.population"],
type = "l", col = "red", lwd = 2,
xlab = "Time (days)", ylab = "Recovering Population (people)",
main = "Recovering Population: 0% vs 100% Antibiotics",
ylim = c(0, max(out_with_antibiotics[, "recovering.population"])))
lines(out_with_antibiotics[, "time"], out_with_antibiotics[, "recovering.population"],
col = "blue", lwd = 2)
legend("bottomright", legend = c("0% Antibiotics", "100% Antibiotics"),
col = c("red", "blue"), lwd = 2)
grid()
# 4. Comparación de DECEASED POPULATION
plot(out_no_antibiotics[, "time"], out_no_antibiotics[, "deceased.population"],
type = "l", col = "red", lwd = 2,
xlab = "Time (days)", ylab = "Deceased Population (people)",
main = "Deceased Population: 0% vs 100% Antibiotics",
ylim = c(0, max(out_no_antibiotics[, "deceased.population"])))
lines(out_with_antibiotics[, "time"], out_with_antibiotics[, "deceased.population"],
col = "blue", lwd = 2)
legend("bottomright", legend = c("0% Antibiotics", "100% Antibiotics"),
col = c("red", "blue"), lwd = 2)
grid()
par(mfrow = c(1, 1))
# Comentarios comparativos:
# INFECTIONS (Nuevas infecciones por día):
# Los picos de infecciones son idénticos (aproximadamente 11 mil personas/día) en ambos escenarios.Los antibióticos no afectan la velocidad de contagio, solo las consecuencias de la infección. Ambas curvas muestran el mismo patrón explosivo inicial.
# DEATHS (Muertes por día):
# Aquí se ve el mayor impacto de los antibióticos. El pico de muertes se reduce de aproximadamente 3 mil a 500 personas/día (reducción del 83%). Las curvas tienen formas similares pero magnitudes completamente diferentes.
# RECOVERING POPULATION (Población recuperada acumulada):
# Los antibióticos invierten completamente el resultado. Sin tratamiento, unas mil personas se recuperan. Con antibióticos, cerca de unas 8,500 personas sobreviven. Es un cambio de 10% a 85% de supervivencia.
# DECEASED POPULATION (Población fallecida acumulada):
# Este es un caso muy crítico: sin antibióticos mueren unas 9 mil personas (90% de la población) versus unas 1,500 personas (15%) con tratamiento. Los antibióticos ayudan a tener un brote manejable versus el escenario de una catastrofe.