HIV Simulation Based on Blower et al. (2000)

Task Overview

This notebook implements a compartmental model of HIV transmission and treatment, adapted from Blower et al. (2000), to simulate and evaluate the impact of antiretroviral therapy (ART) under different assumptions.

Objectives:

• Develop a mathematical model of HIV transmission including ART, drug resistance, and treatment dropout.
• Simulate three scenarios over 20 years:
  • No ART
  • ART (Optimistic): stable behavior, slow resistance emergence
  • ART (Pessimistic): increasing risky behavior, rapid resistance emergence
• Visualize the effects of ART on cumulative HIV infections and deaths.
• Interpret and present results in tabular and graphical formats.

Citation:
Blower, S. M., Gershengorn, H. B., & Grant, R. M. (2000). A tale of two futures: HIV and antiretroviral therapy in San Francisco. Science, 287(5453), 650–654.
https://doi.org/10.1126/science.287.5453.650

---

Summary of the Paper

Blower et al. (2000) use a mathematical model to estimate the population-level impacts of ART on HIV transmission and AIDS-related deaths in the San Francisco gay community. At the time, ~30% of the population was HIV-infected, and ~50% of those infected were receiving ART.

The model incorporates:

• Drug-sensitive and drug-resistant HIV strains
• ART effects on survival and infectivity
• Behavior-driven risk (e.g., changes in partner acquisition rate)
• Resistance emergence and reversion dynamics

Two uncertainty-based scenarios were modeled:

• Optimistic: Low resistance, stable behavior, high ART efficacy
• Pessimistic: High resistance emergence, increased risky behavior

Key results showed that while ART reduces deaths and can lower incidence, increases in risky behavior or drug resistance can significantly reduce its effectiveness. Effective ART rollout must be paired with behavioral interventions.

Model Setup Explanation

This notebook implements an extended SIR-like compartmental model. It tracks five groups over time:

• X: Susceptible individuals
• YSU: Infected with drug-sensitive strain, untreated
• YST: Infected with drug-sensitive strain, treated
• YRU: Infected with drug-resistant strain, untreated
• YRT: Infected with drug-resistant strain, treated
• D: Cumulative AIDS-related deaths (added as an explicit state)

Each individual may:

• Become infected based on contact with infectious individuals
• Start or drop ART
• Develop drug resistance
• Revert from drug-resistant to drug-sensitive status
• Progress to AIDS and die (with progression rate affected by ART)

Model Equations and Definitions

Total Population

The total population is:

$$ N = X + Y_{SU} + Y_{ST} + Y_{RU} + Y_{RT} $$

Forces of Infection

Let $\lambda_S$ and $\lambda_R$ denote the force of infection from drug-sensitive and drug-resistant strains, respectively:

$$ \lambda_S = \frac{\beta_{SU} Y_{SU} + \beta_{ST} Y_{ST}}{N} \quad \quad \lambda_R = \frac{\beta_{RU} Y_{RU} + \beta_{RT} Y_{RT}}{N} $$

Differential Equations

The system is governed by the following ordinary differential equations:

\begin{aligned} \frac{dX}{dt} &= pai - X \left[ c(\lambda_S + \lambda_R) + \mu \right] \\ \frac{dY_{SU}}{dt} &= X c \lambda_S + Y_{RU} q + Y_{ST} g_S - Y_{SU}(s_S + v_{SU} + \mu) \\ \frac{dY_{ST}}{dt} &= Y_{SU} s_S - Y_{ST}(g_S + r + v_{ST} + \mu) \\ \frac{dY_{RU}}{dt} &= X c \lambda_R + Y_{RT} g_R - Y_{RU}(q + e s_R + v_{RU} + \mu) \\ \frac{dY_{RT}}{dt} &= Y_{RU} e s_R + Y_{ST} r - Y_{RT}(g_R + v_{RT} + \mu) \end{aligned}

Parameter Notation

Symbol	Meaning
$pai$	Entry rate into sexually active population
$c$	Number of new partners per year
$\mu$	Natural exit rate
$\beta_{ij}$	Transmission rate for group $i$ (S or R), status $j$ (U or T)
$v_{ij}$	Progression rate for group $i$, status $j$
$s_S$, $s_R$	ART initiation rate for sensitive and resistant strains
$g_S$, $g_R$	ART dropout rates
$r$	Resistance emergence rate
$q$	Reversion rate from drug-resistant to drug-sensitive
$e$	ART efficacy in treating resistant infections

---

import numpy as np
import pandas as pd
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import seaborn as sns

# --- Model Function ---
def blower_model(y, t, params):
    X, YSU, YST, YRU, YRT, D = y
    N = X + YSU + YST + YRU + YRT

    bSU, bST, bRU, bRT = params["bSU"], params["bST"], params["bRU"], params["bRT"]
    vSU, vST, vRU, vRT = params["vSU"], params["vST"], params["vRU"], params["vRT"]
    sS, sR, gS, gR = params["sS"], params["sR"], params["gS"], params["gR"]
    r, q, e, c, pai, mu = params["r"], params["q"], params["e"], params["c"], params["pai"], params["mu"]

    lambdaS = (bSU * YSU + bST * YST) / N
    lambdaR = (bRU * YRU + bRT * YRT) / N

    dX = pai - X * (c * (lambdaS + lambdaR) + mu)
    dYSU = X * c * lambdaS + YRU * q + YST * gS - YSU * (sS + vSU + mu)
    dYST = YSU * sS - YST * (gS + r + vST + mu)
    dYRU = X * c * lambdaR + YRT * gR - YRU * (q + e * sR + vRU + mu)
    dYRT = YRU * e * sR + YST * r - YRT * (gR + vRT + mu)
    dD = YSU * vSU + YST * vST + YRU * vRU + YRT * vRT

    return [dX, dYSU, dYST, dYRU, dYRT, dD]

Stochastic Scenario Simulation

To evaluate the robustness of ART under varying behavioral and clinical conditions, we conducted 100 simulations for each of three policy-relevant scenarios:

• No ART: Baseline scenario with no treatment access
• ART – Optimistic: Low resistance emergence, low dropout, and stable partner behavior
• ART – Pessimistic: High resistance emergence, high dropout, and increased partner acquisition

Randomized Parameter Sampling

For each scenario, key parameters were drawn randomly from plausible ranges to reflect uncertainty in:

• Behavioral factors (e.g., partner acquisition rate c, dropout rate gS, gR)
• Clinical dynamics (e.g., resistance emergence r, ART effectiveness bST)
• Programmatic features (e.g., ART initiation sS, sR)

This approach creates probabilistic envelopes around model outcomes rather than relying on fixed-point estimates.

Outcome Interpretation

We plotted the median, 25th percentile, and 75th percentile for total HIV infections over 20 years across scenarios. These uncertainty bands help illustrate:

• The effectiveness and variability of ART under real-world constraints
• The potential for backfire effects in pessimistic settings if resistance and risky behavior are not controlled
• The comparative public health gains achievable under favorable conditions

# --- Initial Conditions & Time ---
init = [7000, 1000, 500, 300, 200, 0]
t = np.linspace(0, 20, 201)  # 20 years, 0.1-year steps

# --- Randomized Parameter Generator ---
def generate_params_for_scenario(scenario):
    base = {
        "pai": 2133,
        "bSU": 0.1,
        "bRU": 0.05,
        "bRT": 0.01,
        "vSU": 1/12,
        "vST": 1/24,
        "vRU": 1/14,
        "vRT": 1/18,
        "q": 1/0.15,
        "e": 0.8,
        "mu": 1/30,
    }

    if scenario == "No ART":
        return {
            **base,
            "bST": 0.01,
            "c": np.random.uniform(1.5, 2.5),
            "sS": 0,
            "sR": 0,
            "r": 0,
            "gS": 0,
            "gR": 0,
        }
    elif scenario == "ART - Optimistic":
        return {
            **base,
            "bST": np.random.uniform(0.005, 0.015),
            "c": np.random.uniform(1.5, 1.8),
            "sS": 0.1,
            "sR": 0.05,
            "r": np.random.uniform(0.05, 0.15),
            "gS": np.random.uniform(0.03, 0.07),
            "gR": np.random.uniform(0.08, 0.12),
        }
    elif scenario == "ART - Pessimistic":
        return {
            **base,
            "bST": np.random.uniform(0.03, 0.06),
            "c": np.random.uniform(2.0, 2.7),
            "sS": 0.1,
            "sR": 0.05,
            "r": np.random.uniform(0.4, 0.6),
            "gS": np.random.uniform(0.12, 0.18),
            "gR": np.random.uniform(0.22, 0.28),
        }

# --- Simulate with and without ART ---
def simulate_randomized_scenarios(scenarios, num_runs=100):
    all_data = []
    for scenario in scenarios:
        for run_id in range(num_runs):
            params = generate_params_for_scenario(scenario)
            sol = odeint(blower_model, init, t, args=(params,))
            df = pd.DataFrame(sol, columns=["X", "YSU", "YST", "YRU", "YRT", "D"])
            df["TotalInfected"] = df[["YSU", "YST", "YRU", "YRT"]].sum(axis=1)
            df["Scenario"] = scenario
            df["time"] = t
            df["run"] = run_id
            all_data.append(df)
    return pd.concat(all_data, ignore_index=True)

# --- Run Simulations ---
scenarios = ["No ART", "ART - Optimistic", "ART - Pessimistic"]
results = simulate_randomized_scenarios(scenarios, num_runs=100)

# --- Compute uncertainty bands ---
summary = (
    results.groupby(["Scenario", "time"])["TotalInfected"]
    .quantile([0.25, 0.5, 0.75])
    .unstack()
    .rename(columns={0.25: "q25", 0.5: "q50", 0.75: "q75"})
    .reset_index()
)

# --- Plot Uncertainty Band for Total Infected ---
plt.figure(figsize=(10, 6))
colors = {"No ART": "red", "ART - Optimistic": "green", "ART - Pessimistic": "green"}
linestyles = {"No ART": "-", "ART - Optimistic": "-", "ART - Pessimistic": "--"}

for scenario in scenarios:
    df = summary[summary["Scenario"] == scenario]
    plt.plot(df["time"], df["q50"], label=f"{scenario} (median)", color=colors[scenario], linestyle=linestyles[scenario])
    plt.fill_between(df["time"], df["q25"], df["q75"], alpha=0.2, color=colors[scenario])

plt.xlabel("Time (years)")
plt.ylabel("Total Infected")
plt.title("HIV Simulation with Randomized Scenario Parameters")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

d_summary = (
    results.groupby(["Scenario", "time"])["D"]
    .quantile([0.25, 0.5, 0.75])
    .unstack()
    .rename(columns={0.25: "q25", 0.5: "q50", 0.75: "q75"})
    .reset_index()
)

# --- Plot Uncertainty Band for Total Death ---
plt.figure(figsize=(10, 6))
colors = {"No ART": "red", "ART - Optimistic": "green", "ART - Pessimistic": "green"}
linestyles = {"No ART": "-", "ART - Optimistic": "-", "ART - Pessimistic": "--"}

for scenario in scenarios:
    df = d_summary[d_summary["Scenario"] == scenario]
    plt.plot(df["time"], df["q50"], label=f"{scenario} (median)", color=colors[scenario], linestyle=linestyles[scenario])
    plt.fill_between(df["time"], df["q25"], df["q75"], alpha=0.2, color=colors[scenario])

plt.xlabel("Time (years)")
plt.ylabel("Total Deaths")
plt.title("HIV Simulation with Randomized Scenario Parameters")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Summary of HIV Simulation Results

Total Infections Over 20 Years

• ART (Optimistic):
  Case counts remain relatively low (~1736), indicating that ART with low dropout and stable behavior can contain the epidemic effectively.
  

• No ART:
  Infections rise steadily to ~4523, showing clear but slower epidemic expansion in the absence of intervention.
  

• ART (Pessimistic):
  Infections escalate to over 6000, suggesting that ART without behavioral control and with high resistance can worsen the epidemic.

• Insight: ART policies must incorporate behavioral and resistance safeguards to prevent counterproductive outcomes.

---

Cumulative Deaths Over 20 Years

• ART (Optimistic):
  Deaths are held to ~2510, showcasing ART’s ability to prolong life when implemented well.

• No ART:
  Deaths reach ~4542, representing the burden of untreated HIV progression.

• ART (Pessimistic):
  Deaths climb to ~5441, worse than no ART, driven by resistance, dropout, and increased transmission.

• Insight: ART without effective adherence and behavioral strategies can lead to higher mortality than no intervention.

---

Overall Summary Table

Scenario	Infections (Year 20)	Deaths (Year 20)	Key Takeaway
ART - Optimistic	~1736	~2510	Best-case outcome, stable epidemic
No ART	~4523	~4542	Moderate growth and mortality
ART - Pessimistic	~6005	~5441	Worst-case; ART without safeguards can worsen HIV

Alignment with Blower et al. (2000)

This simulation reproduces the core findings of Blower et al. (2000):
ART can dramatically reduce infections and deaths—but only if resistance emergence and risky behavior are minimized.
Otherwise, it may amplify rather than suppress the epidemic.