{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c85766b3-356b-46a9-a6cf-931ba307f1c1",
   "metadata": {},
   "source": [
    "# Vector-borne disease"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e11e1d1-8535-44b8-867b-9a7f338d01ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pykappa import System"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd1c8f6d-07e9-4dee-8ec0-a8df7fa84eec",
   "metadata": {},
   "source": [
    "The following model is adapted from [\"Rule-based epidemic models\"](https://doi.org/10.1016/j.jtbi.2021.110851) by Waites et al. We simulate a simple vector-borne disease which can travel between both mosquitoes (`V(x{s i})`) and hosts (`P(x{s e i r})`). The disease progresses: hosts are initially susceptible (`x{s}`) to the disease, later become exposed after coming in contact with a mosquito (`x{e}`), are infected (`x{i}`), and eventually recover (`x{r}`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc1a16de-fc90-4ac9-a93a-1858d7d81543",
   "metadata": {},
   "outputs": [],
   "source": [
    "INIT_M=50_000           # intial number of mosquitoes;\n",
    "INIT_N=10_000           # total host population, of which...\n",
    "INIT_I=100              # are initially infectious, and\n",
    "INIT_S=INIT_N - INIT_I  # are initially susceptible\n",
    "\n",
    "variables = {\n",
    "    \"beta\": \"0.036\",    # probability of infection from a bite\n",
    "    \"bprime\": \"1\",      # proabability of a vector becoming infectious\n",
    "    \"alpha\": \"0.2\",     # progression from exposed to infectious\n",
    "    \"gamma\": \"0.1429\",  # progression from infectious to removed\n",
    "    \"kappa\": \"1.0\",     # bites per day per mosquito\n",
    "    \"kb\": \"0.1429\",     # birth rate\n",
    "    \"kd\": \"0.1429\",     # death rate\n",
    "    \"water\": \"1\"        # amount of breeding habitat available\n",
    "}\n",
    "\n",
    "rules = [\n",
    "    \"P(x{e}) -> P(x{i}) @ 'alpha'\",  # disease progression\n",
    "    \"P(x{i}) -> P(x{r}) @ 'gamma'\",  # disease recovery\n",
    "    \n",
    "    \"V(), . -> V(), V(x{s}) @ 'water' * 'kb'\",  # vector birth\n",
    "    \"V()    -> .            @ 'kd'\",            # vector annihilation\n",
    "\n",
    "    # both hosts and vectors can be infected from a bite event\n",
    "    \"P(x{i}), V(x{s}) -> P(x{i}), V(x{i}) @ 'bprime' * 'kappa' / 'M'\",\n",
    "    \"V(x{i}), P(x{s}) -> V(x{i}), P(x{e}) @ 'beta' * 'kappa' / 'N'\"\n",
    "]\n",
    "\n",
    "observables = {\n",
    "    \"N\": \"|P()|\",\n",
    "    \"S\": \"|P(x{s})|\",\n",
    "    \"E\": \"|P(x{e})|\",\n",
    "    \"I\": \"|P(x{i})|\",\n",
    "    \"R\": \"|P(x{r})|\",\n",
    "    \n",
    "    \"M\": \"|V()|\",\n",
    "    \"Vs\": \"|V(x{s})|\",\n",
    "    \"Vi\": \"|V(x{i})|\"\n",
    "}\n",
    "\n",
    "system = System.from_kappa(\n",
    "    mixture={\"P(x{i})\": INIT_I, \"P(x{s})\": INIT_S, \"V(x{s})\": INIT_M},\n",
    "    rules=rules,\n",
    "    variables=variables,\n",
    "    observables=observables,\n",
    "    seed=42\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef4bbb3a-f51c-4960-8e57-d047e82ab931",
   "metadata": {},
   "source": [
    "We run the simulation for 375 time units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4fa3792-9590-44d4-8b57-4130449aad16",
   "metadata": {},
   "outputs": [],
   "source": [
    "while system.time < 375:\n",
    "    system.update()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98032bac-b30a-4801-a06e-ee1ceeb5772e",
   "metadata": {},
   "source": [
    "We plot the numbers of hosts each disease progression state over time. The system eventually equilibriates; the disease is wiped out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2090b327-db7c-4fae-859d-6d507cb4fd5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.monitor.plot(combined=True, observables=[\"S\", \"E\", \"I\", \"R\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfe32cc1-8db4-4487-b6c1-ef6daf262341",
   "metadata": {},
   "source": [
    "All infected vectors eventually perish:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a5e937d-0535-4508-a748-9f0e1823d7d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.monitor.plot(combined=True, observables=[\"Vs\", \"Vi\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74d6108f-b298-467c-a0c4-889e0459e6b8",
   "metadata": {},
   "source": [
    "Now we regenerate the system from scratch in order to demonstrate the impact of an environmental intervention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ebda3f3-6cc3-47c6-b091-707fc463affc",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = System.from_kappa(\n",
    "    mixture={\"P(x{i})\": INIT_I, \"P(x{s})\": INIT_S, \"V(x{s})\": INIT_M},\n",
    "    rules=rules,\n",
    "    variables=variables,\n",
    "    observables=observables,\n",
    "    seed=42\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ea0aa04-8bb1-454c-8913-fbfd1b8ea49a",
   "metadata": {},
   "source": [
    "Every sixty time units, we remove water from the environment, which slows the mosquito birth rate relative to the death rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3deb2ef0-b63f-40fd-babf-51f5fa30626e",
   "metadata": {},
   "outputs": [],
   "source": [
    "update_after = 60.0\n",
    "while system[\"M\"] > 0:\n",
    "    system.update()\n",
    "\n",
    "    if system.time > update_after:\n",
    "        update_after += 60.0\n",
    "        system[\"water\"] = system[\"water\"] * 0.9"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a514ec6-7a07-41ff-8db7-f1f0197aeddb",
   "metadata": {},
   "source": [
    "All mosquitoes eventually perish, and the disease is eliminated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24439af6-f9ff-4a8d-8558-00312e732248",
   "metadata": {},
   "outputs": [],
   "source": [
    "system.monitor.plot(combined=True, observables=[\"S\", \"E\", \"I\", \"R\"]);\n",
    "system.monitor.plot(combined=True, observables=[\"Vs\", \"Vi\"]);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}