30:T598,md`# Epidemic Simulation This notebook implements a [simple infectious disease model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). There are ${N} inviduals who are either Susceptible ${box(STATES.SUSCEPTIBLE)}, Infected ${box(STATES.INFECTED)}, or Recovered ${box(STATES.RECOVERED)}. They are connected to each other in a network generated using the [Barabási–Albert preferential attachment model](https://observablehq.com/@herbps10/barabasi-albert-model). The disease spreads between individuals that are connected in the network. The simulation works like this: 1. ${num_infected_initialization} individuals are chosen at random to be infected at the beginning of the simulation. All other indivduals are susceptible. 2. Infected individuals have a ${Math.round(infection_probability * 100)}% probability of infecting each of their susceptible neighbors at each timestep. 3. After ${time_infected} timesteps, infected inviduals move to the recovered state. 4. After ${time_recovered} timesteps, recovered individuals return to the susceptible state. A few things to try: - *Single epidemic*: If the probability of infection is high, then the epidemic can spread quickly across the network. The epidemic will quickly use up the susceptible individuals and die out. - *Recurring epidemics*: If the recovery time is less than the infected time, you can induce an epidemic that continues indefinitely. `31:Tab4,chart = { const xExtent = d3.extent(graph.nodes.map((d) => d.x)); const yExtent = d3.extent(graph.nodes.map((d) => d.y)); elements.svg.attr("viewBox", `${xExtent[0] - 5} ${yExtent[0] - 5} ${xExtent[1] - xExtent[0] + 5} ${yExtent[1] - yExtent[0] + 10}`); elements.link .selectAll("line") .data(graph.links) .join("line") .attr("stroke-width", d => Math.sqrt(d.value)) .attr("x1", d => d.source.x) .attr("y1", d => d.source.y) .attr("x2", d => d.target.x) .attr("y2", d => d.target.y); elements.link_spreading .selectAll("line") .data(graph.links) .join("line") .attr("stroke-width", d => Math.sqrt(d.value)) .attr("stroke", color(STATES.INFECTED)) .attr("stroke-opacity", 0.6) .attr("stroke-dasharray", 1000) .attr("x1", d => d.source.x) .attr("y1", d => d.source.y) .attr("x2", d => d.target.x) .attr("y2", d => d.target.y) const spread_direction = (d) => { const direction1 = simulation_results[visualization_timestep].infection_pathways.filter((pathway) => (pathway.source == d.source.index && pathway.target == d.target.index)); const direction2 = simulation_results[visualization_timestep].infection_pathways.filter((pathway) => (pathway.source == d.target.index && pathway.target == d.source.index)); return (direction1.length > 0) ? 1 : (direction2.length > 0) ? -1 : 0; } elements.link_spreading .selectAll("line") .attr("stroke-dashoffset", (d) => { if(spread_direction(d) == 1) return 1000; if(spread_direction(d) == -1) return -1000; return 1000; }) .transition() .duration(speed / 2) .attr("stroke-dashoffset", (d) => { // Check if the infection spread across this edge const matches = simulation_results[visualization_timestep].infection_pathways.filter((pathway) => (pathway.source == d.source.index && pathway.target == d.target.index) || (pathway.source == d.target.index && pathway.target == d.source.index) ); return matches.length > 0 ? 0 : 1000; }) .transition() .duration(speed / 2) .attr("stroke-dashoffset", (d) => { if(spread_direction(d) == 1) return -1000; if(spread_direction(d) == -1) return 1000; return 1000; }); elements.node .selectAll("circle") .data(graph.nodes) .join("circle") .attr("r", 4) .attr("cx", d => d.x) .attr("cy", d => d.y) .transition() .duration(250) .delay(speed / 2 - 100) .attr("fill", (d, i) => color(simulation_results[visualization_timestep].node_states[i])); invalidation.then(() => simulation.stop()); return elements.svg.node(); }32:T71f,{ const margin = { top: 20, right: 30, bottom: 30, left: 40 }; const height = width / 3; const svg = d3.create("svg") .attr("viewBox", [0, 0, width, height]); const xAxis = g => g .attr("transform", `translate(0,${height - margin.bottom})`) .call(d3.axisBottom(x).ticks(width / 80).tickSizeOuter(0)) const yAxis = g => g .attr("transform", `translate(${margin.left},0)`) .call(d3.axisLeft(y)) .call(g => g.select(".domain").remove()) .call(g => g.select(".tick:last-of-type text").clone() .attr("x", 0) .attr("y", -15) .attr("text-anchor", "start") .attr("font-weight", "bold") .text("Number of individuals")) const x = d3.scaleLinear() .domain([0, timesteps]) .range([margin.left, width - margin.right]); const y = d3.scaleLinear() .domain([0, N]) .range([height - margin.bottom, margin.top]) const line = d3.line() .defined(d => !isNaN(d)) .x((d, i) => x(i)) .y(d => y(d)) svg.append("g") .call(xAxis); svg.append("g") .call(yAxis); svg.selectAll("path.timeseries") .data(summary_statistics.state_timeseries) .enter() .append("path") .attr("class", "timeseries") .attr("fill", "none") .attr("stroke", (d) => color(d.state)) .attr("stroke-width", 1.5) .attr("stroke-linejoin", "round") .attr("stroke-linecap", "round") .attr("d", (d) => line(d.values)); svg.selectAll("line.scrubber") .data([visualization_timestep]) .enter() .append("line") .attr("class", "scrubber") .attr("x1", (d) => x(d)) .attr("x2", (d) => x(d)) .attr("y1", margin.top) .attr("y2", height - margin.bottom) .attr("stroke", "#ccc") .attr("stroke-dasharray", "2 2"); return svg.node(); }33:T5c4,tick_simulation = function(state, edge_list) { const new_state = { node_states: state.node_states.slice(0), time_in_state: state.time_in_state.slice(0), infection_pathways: []}; // Update time spent in state for(let i = 0; i < edge_list.length; i++) { new_state.time_in_state[i] = state.time_in_state[i] + 1; } // Loop through each individual for(let i = 0; i < edge_list.length; i++) { if(state.node_states[i] == STATES.INFECTED) { // If they have been infected for N timesteps, set them back to susceptible if(state.time_in_state[i] == time_infected) { new_state.node_states[i] = STATES.RECOVERED; new_state.time_in_state[i] = 0; } // If they are still infected, infect their neighbors for(let j = 0; j < edge_list[i].length; j++) { if(state.node_states[edge_list[i][j]] == STATES.SUSCEPTIBLE && Math.random() < infection_probability) { new_state.node_states[edge_list[i][j]] = STATES.INFECTED; new_state.time_in_state[edge_list[i][j]] = 0; new_state.infection_pathways.push({ source: i, target: edge_list[i][j] }); } } } if(state.node_states[i] == STATES.RECOVERED && state.time_in_state[i] == time_recovered) { new_state.node_states[i] = STATES.SUSCEPTIBLE; new_state.time_in_state[i] = 0; } } return new_state; }34:[] 2e:{"title":"Epidemic Simulation","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":390,"value":"viewof rerun_simulation = button({\n value: \"Generate new epidemic\",\n description: \"Re-initialize the simulation.\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":271,"value":"viewof visualization_timestep = Scrubber(Array.from({length: timesteps}, (_, i) => i), {\n delay: speed + 50,\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"viewof time_infected = slider({\n min: 1,\n max: 10,\n value: 3,\n step: 1,\n title: \"Time infected\",\n description: \"How many timesteps an individual stays in the infected state.\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":372,"value":"viewof time_recovered = slider({\n min: 1,\n max: 10,\n value: 2,\n step: 1,\n title: \"Time recovered\",\n description: \"How many timesteps an individual stays in the recovered state.\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":386,"value":"viewof infection_probability = slider({\n min: 0,\n max: 1,\n title: \"Probability of spreading infection\",\n description: \"Probability that an infected individual infects a susceptible neighbor in one timestep.\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":532,"value":"md`### Epidemic dynamics`","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":529,"value":"md`## More settings`","pinned":false,"hidden":false,"mode":"ojs"},{"id":110,"value":"viewof N = slider({\n title: \"Number of individuals\",\n min: 1,\n max: 500,\n step: 1,\n value: 200,\n description: \"Size of the network.\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":80,"value":"viewof attachment_strength = slider({\n title: \"Preferential attachment strength\",\n min: 0.1,\n max: 2,\n value: 1,\n description: \"Controls the topology of the network. Higher attachment strength induces a graph where a few individuals have many neighbors, and most individuals only have one or two neighbors.\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":610,"value":"viewof num_infected_initialization = slider({\n min: 1,\n max: Math.min(N, 5),\n step: 1,\n title: \"Number of individuals to infect at initialization\",\n value: 2\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":251,"value":"viewof timesteps = slider({\n min: 10,\n max: 100,\n step: 10,\n title: \"Number of timesteps\",\n value: 50\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":88,"value":"viewof rerun = button({ value: \"Generate new network\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":592,"value":"md`\nPossible extensions:\n- *Implement different network structures*, like a small-world network or regular lattice.\n- *Vaccinations*: allow individuals to be vaccinated, protecting them from infection.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":694,"value":"speed = 1000","pinned":false,"hidden":false,"mode":"ojs"},{"id":555,"value":"md`## Simulation code`","pinned":true,"hidden":false,"mode":"ojs"},{"id":241,"value":"simulation_results = {\n rerun_simulation;\n const states = [];\n states[0] = initialize_simulation(edge_list(graph));\n for(let timestep = 1; timestep < timesteps; timestep++) {\n states[timestep] = tick_simulation(states[timestep - 1], edge_list(graph));\n }\n \n return states;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":255,"value":"initialize_simulation = function(edge_list) {\n const state = {\n node_states: [],\n time_in_state: [],\n infection_pathways: []\n };\n \n // Set all nodes to susceptible\n for(let i = 0; i < edge_list.length; i++) {\n state.node_states[i] = STATES.SUSCEPTIBLE; \n state.time_in_state[i] = 0;\n }\n \n // Pick random nodes to be infected\n const shuffled = Array.from({length: N}, (_, i) => i).sort(() => 0.5 - Math.random());\n const selected = shuffled.slice(0, num_infected_initialization);\n for(let index in selected) {\n state.node_states[index] = STATES.INFECTED;\n state.time_in_state[index] = 0;\n }\n return state;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":257,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":397,"value":"summary_statistics = {\n const num_in_states = [];\n for(let i = 0; i < timesteps; i++) {\n num_in_states[i] = count_states(simulation_results[i]);\n }\n \n return ({\n state_timeseries: [\n { name: \"Susceptible\", state: STATES.SUSCEPTIBLE, values: num_in_states.map((d) => d[STATES.SUSCEPTIBLE]) },\n { name: \"Infected\", state: STATES.INFECTED, values: num_in_states.map((d) => d[STATES.INFECTED]) },\n { name: \"Recovered\", state: STATES.RECOVERED, values: num_in_states.map((d) => d[STATES.RECOVERED]) }\n ]\n });\n};","pinned":false,"hidden":false,"mode":"ojs"},{"id":398,"value":"count_states = (state) => {\n const states = [0, 0, 0];\n for(let i = 0; i < state.node_states.length; i++) {\n states[state.node_states[i]] += 1;\n }\n \n return states;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":382,"value":"barabasi_num_links = 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"elements = {\n const svg = d3.select(DOM.svg(width, height));\n const link = svg.append(\"g\")\n .attr(\"stroke\", \"#999\")\n .attr(\"stroke-opacity\", 0.6);\n \n const link_spreading = svg.append(\"g\");\n \n const node = svg.append(\"g\")\n .attr(\"stroke\", \"#fff\")\n .attr(\"stroke-width\", 1.5);\n \n return ({\n svg,\n link,\n link_spreading,\n node\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":283,"value":"simulation = d3.forceSimulation(graph.nodes)\n .force(\"link\", d3.forceLink(graph.links).id(d => d.id))\n .force(\"charge\", d3.forceManyBody().strength(-6))\n .force(\"center\", d3.forceCenter(width / 2, height / 2))\n .tick(1000)","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"graph = {\n rerun;\n return barabasi_albert(N, attachment_strength, barabasi_num_links)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":242,"value":"edge_list = (graph) => {\n const list = [];\n for(let i = 0; i < graph.links.length; i++) {\n if(list[graph.links[i].source.id] == null) list[graph.links[i].source.id] = [];\n if(list[graph.links[i].target.id] == null) list[graph.links[i].target.id] = [];\n \n list[graph.links[i].source.id].push(graph.links[i].target.id);\n list[graph.links[i].target.id].push(graph.links[i].source.id);\n }\n \n return list;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":464,"value":"box = function(state) {\n return html`\n \n `;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":261,"value":"STATES = ({\n SUSCEPTIBLE: 0,\n INFECTED: 1,\n RECOVERED: 2\n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":18,"value":"height = width * 0.5","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"barabasi_albert = function(n, attachment_strength, num_links) {\n const nodes = [{ id: 0, degree: 0 }];\n const links = [];\n for(let i = 1; i < n; i++) {\n // Preferentially attach node to existing nodes\n // with high degree\n const indices = d3.range(0, i);\n \n // We weight which node to attach to by the node's degree\n const weights = nodes\n .map((d) => nodes.length == 1 ? 1 : d.degree)\n .map((d) => Math.pow(d, attachment_strength));\n \n for(let j = 0; j < num_links; j++) {\n const attachment_node = weighted_random(indices, weights);\n\n // Add the link\n links.push({ source: i, target: attachment_node });\n nodes[attachment_node].degree += 1;\n }\n \n // Add the new node\n nodes.push({ id: i, degree: 1 });\n }\n \n return ({ nodes, links });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":33,"value":"weighted_random = function(x, weights) {\n const weight_sum = d3.sum(weights);\n const normalized_weights = weights.map((d) => d / weight_sum);\n \n let r = Math.random();\n for(let i = 0; i < x.length; i++) {\n if(r < normalized_weights[i]) {\n return x[i];\n }\n else {\n r -= normalized_weights[i];\n }\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"drag = simulation => {\n \n function dragstarted(d) {\n if (!d3.event.active) simulation.alphaTarget(0.3).restart();\n d.fx = d.x;\n d.fy = d.y;\n }\n \n function dragged(d) {\n d.fx = d3.event.x;\n d.fy = d3.event.y;\n }\n \n function dragended(d) {\n if (!d3.event.active) simulation.alphaTarget(0);\n d.fx = null;\n d.fy = null;\n }\n \n return d3.drag()\n .on(\"start\", dragstarted)\n .on(\"drag\", dragged)\n .on(\"end\", dragended);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":274,"value":"color = (node_state) => {\n if(node_state == STATES.INFECTED) return palette.pomegranate;\n if(node_state == STATES.RECOVERED) return palette.belize_hole;\n if(node_state == STATES.SUSCEPTIBLE) return palette.concrete;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":21,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":511,"value":"import {Scrubber} from \"@mbostock/scrubber\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":543,"value":"import { palette } from \"@herbps10/flat-ui-colors\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":76,"value":"import { slider, button } from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}