Epidemics | Computational Physics

Introduction

This is a very brief post, with a simulation written in javascript. It’s a very simple model of how a dissease can spread into a population. You may get the code, adjust the parameters and watch different outcomes.
The idea of this project came from this: Why outbreaks like coronavirus spread exponentially, and how to ‘flatten the curve’. I highly recommend reading it.

Theory

I recommend watching this video, for some theory:

It’s a good start, there is a lot more information about this on the internet, so I won’t bother with formulae here.

Simulation

The model resembles a lot the one from the link pointed above, with an addition: this model also simulates death, not all ‘people’ recover.

Briefly, it’s just molecular dynamics in 2D, in a ‘box’, with elastic collisions. All ‘particles’ have the same mass and they can be either not infected, infected, cured and dead. The dead ones disappear from the simulation.
Once an infected ‘particle’ hits one that wasn’t infected yet, it will infect it. A ‘particle’ once having the illness, has a chance of dying (but only in the last two thirds of the infected period of time). After being sick for a while, if it doesn’t die in the meantime, it will recover and have immunity, so it cannot be infected again.

There is also a chart showing the number of infections over time. You’ll notice that it will level out under 100%, the difference being the ones that did not get an infection, being lucky, and the dead ones.

The code has a feature that it’s usually a bug: because of the too big time step, some ‘particles’ can stick together. Basically when the collide they overlap a little (more, if the time step is big) and if the cannot go apart in one time step, remaining overlapped, they will be considered as colliding again and so on. There are various ways of solving this but for this model I consider it a feature: people sometimes do something analogous.

The blue ‘particles’ are not infected, the red ones are infected and the green ones are cured. The ones that die disappear from the canvas.

0

The computation will be a little slow on a phone, for this you should better use a desktop.

Related to this, for a better implementation on such molecular dynamics where there is short range interaction only, visit Event Driven Molecular Dynamics. It’s also 3D and with nicer graphics, using OpenGL.

The code

There are some comments in the code, hopefully it’s clear enough. I wrote it very fast, so it’s far from perfect but it seems to work as intended.

(function () {
    var canvas = document.getElementById("Canvas");
    var chart = document.getElementById("Chart");

    canvas.width = Math.min(window.innerWidth - 50, window.innerHeight - 50);
    if (canvas.width < 300) canvas.width = 300;
    else if (canvas.width > 800) canvas.width = 800;
    canvas.height = canvas.width;

    chart.width = canvas.width;
    chart.height = chart.width / 3;

    var ctx = canvas.getContext("2d");
    var canvasText = document.getElementById("canvasText");

    var ctxChart = chart.getContext("2d");
    ctxChart.strokeStyle = "#FF0000";
    ctxChart.lineWidth = 4;

    var people = [];

    // some parameters
    var nrPeople = 500;
    var speed = canvas.width / 50.;
    var radius = canvas.width / 100.;
    var deltat = 0.01;
    var cureTime = 5.;
    var deathProb = 0.05;
    var infectProb = 1.;

    // chart
    var chartPosX = 0;
    var chartValY = 0;
    var chartXScale = 0.1;
    var chartYScale = chart.height / nrPeople;

    // statistics
    var deaths = 0;
    var infections = 1;
    var cures = 0;

    function Init() {
        deaths = 0;
        infections = 1;
        cures = 0;
        people = [];
        chartPosX = 0;
        chartValY = 0;

        for (i = 0; i < nrPeople; ++i) {
            var person =
            {
                posX: 0,
                posY: 0,
                velX: 0,
                velY: 0,
                dead: false,
                infected: false,
                cured: false,
                infectedTime: 0.0,

                NormalizeVelocity: function () {
                    var len = Math.sqrt(this.velX * this.velX + this.velY * this.velY);
                    this.velX /= len;
                    this.velY /= len;
                },

                Distance: function (other) {
                    var distX = this.posX - other.posX;
                    var distY = this.posY - other.posY;

                    return Math.sqrt(distX * distX + distY * distY);
                },

                Collision: function (other) {
                    var dist = this.Distance(other);

                    return dist <= 2. * radius;
                },

                Collide: function (other) {
                    var velXdif = this.velX - other.velX;
                    var velYdif = this.velY - other.velY;

                    var posXdif = this.posX - other.posX;
                    var posYdif = this.posY - other.posY;

                    var dist2 = posXdif * posXdif + posYdif * posYdif;
                    var dotProd = velXdif * posXdif + velYdif * posYdif;
                    dotProd /= dist2;

                    this.velX -= dotProd * posXdif;
                    this.velY -= dotProd * posYdif;

                    other.velX += dotProd * posXdif;
                    other.velY += dotProd * posYdif;
                }
            };

            for (; ;) {
                var X = Math.floor(Math.random() * (canvas.width - 2. * radius)) + radius;
                var Y = Math.floor(Math.random() * (canvas.height - 2 * radius)) + radius;

                person.posX = X;
                person.posY = Y;

                overlap = false;
                for (j = 0; j < i; ++j) {
                    var person2 = people[j];
                    if (person2.Collision(person)) {
                        overlap = true;
                        break;
                    }
                }
                if (!overlap) break;
            }

            person.velX = Math.random() - 0.5;
            person.velY = Math.random() - 0.5;
            person.NormalizeVelocity();
            person.velX *= speed;
            person.velY *= speed;

            if (i == 0) person.infected = true;

            people.push(person);
        }
    }

    Init();

    function Advance() {
        // for each from the population
        for (i = 0; i < nrPeople; ++i) {
            var person = people[i];
            if (person.dead) continue;

            // move

            person.posX += person.velX * deltat;
            person.posY += person.velY * deltat;

            // collide / infect / cure

            // first, walls collisions
            if (person.posX <= radius) {
                person.velX *= -1;
                person.posX = radius;
            }
            else if (person.posX >= canvas.width - radius) {
                person.velX *= -1;
                person.posX = canvas.width - radius;
            }

            if (person.posY <= radius) {
                person.velY *= -1;
                person.posY = radius;
            }
            else if (person.posY >= canvas.height - radius) {
                person.velY *= -1;
                person.posY = canvas.height - radius;
            }

            // keep track of how long the infection lasts
            if (person.infected)
                person.infectedTime += deltat;

            // collisions and infections between them
            for (j = 0; j < i; ++j) {
                var person2 = people[j];
                if (person2.dead) continue;

                if (person.Collision(person2)) {
                    person.Collide(person2);

                    if (person.infected && !person2.infected && !person2.cured) {
                        if (Math.random() < infectProb)
                        {
                           person2.infected = true;
                           ++infections;
                        }
                    }
                    else if (person2.infected && !person.infected && !person.cured) {
                        if (Math.random() < infectProb)
                        {
                           person.infected = true;
                           ++infections;
                        }
                    }
                }
            }

            // cure
            if (person.infected && person.infectedTime > cureTime) {
                person.infected = false;
                person.cured = true;
                ++cures;
            }

            // kill
            if (person.infected && person.infectedTime > cureTime / 3.) {
                if (Math.random() < deathProb * deltat / (cureTime * 2. / 3.)) {
                    person.dead = true;
                    ++deaths;
                }
            }
        }

        // display
        ctx.clearRect(0, 0, canvas.width, canvas.height);
        for (i = 0; i < nrPeople; ++i) {
            var person = people[i];
            if (person.dead) continue;
            ctx.beginPath();
            ctx.arc(person.posX, person.posY, radius, 0, 2 * Math.PI, false);
            ctx.fillStyle = person.infected ? 'red' : (person.cured ? 'green' : 'blue');
            ctx.fill();
            ctx.stroke();
        }

        canvasText.innerHTML = "Total: " + nrPeople + " Infected: " + infections + " Deaths: " + deaths + " Cured: " + cures + " Sick: " + (infections - cures - deaths);

        ctxChart.beginPath();
        ctxChart.moveTo(chartPosX * chartXScale, chart.height - chartValY * chartYScale);
        ++chartPosX;  
        chartValY = infections;      
        ctxChart.lineTo(chartPosX * chartXScale, chart.height - chartValY * chartYScale);
        ctxChart.stroke();
       
        if (infections - deaths == cures) // there is no more left to cure
        {
            ctxChart.clearRect(0, 0, chart.width, chart.height);
            Init();
        }
    }

    setInterval(Advance, 10);
})();

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

(function () {

var canvas = document.getElementById("Canvas");

var chart = document.getElementById("Chart");

canvas.width = Math.min(window.innerWidth - 50, window.innerHeight - 50);

if (canvas.width < 300) canvas.width = 300;

else if (canvas.width > 800) canvas.width = 800;

canvas.height = canvas.width;

chart.width = canvas.width;

chart.height = chart.width / 3;

var ctx = canvas.getContext("2d");

var canvasText = document.getElementById("canvasText");

var ctxChart = chart.getContext("2d");

ctxChart.strokeStyle = "#FF0000";

ctxChart.lineWidth = 4;

var people = [];

// some parameters

var nrPeople = 500;

var speed = canvas.width / 50.;

var radius = canvas.width / 100.;

var deltat = 0.01;

var cureTime = 5.;

var deathProb = 0.05;

var infectProb = 1.;

// chart

var chartPosX = 0;

var chartValY = 0;

var chartXScale = 0.1;

var chartYScale = chart.height / nrPeople;

// statistics

var deaths = 0;

var infections = 1;

var cures = 0;

function Init() {

deaths = 0;

infections = 1;

cures = 0;

people = [];

chartPosX = 0;

chartValY = 0;

for (i = 0; i < nrPeople; ++i) {

var person =

{

posX: 0,

posY: 0,

velX: 0,

velY: 0,

dead: false,

infected: false,

cured: false,

infectedTime: 0.0,

NormalizeVelocity: function () {

var len = Math.sqrt(this.velX * this.velX + this.velY * this.velY);

this.velX /= len;

this.velY /= len;

Distance: function (other) {

var distX = this.posX - other.posX;

var distY = this.posY - other.posY;

return Math.sqrt(distX * distX + distY * distY);

Collision: function (other) {

var dist = this.Distance(other);

return dist <= 2. * radius;

Collide: function (other) {

var velXdif = this.velX - other.velX;

var velYdif = this.velY - other.velY;

var posXdif = this.posX - other.posX;

var posYdif = this.posY - other.posY;

var dist2 = posXdif * posXdif + posYdif * posYdif;

var dotProd = velXdif * posXdif + velYdif * posYdif;

dotProd /= dist2;

this.velX -= dotProd * posXdif;

this.velY -= dotProd * posYdif;

other.velX += dotProd * posXdif;

other.velY += dotProd * posYdif;

}

};

for (; ;) {

var X = Math.floor(Math.random() * (canvas.width - 2. * radius)) + radius;

var Y = Math.floor(Math.random() * (canvas.height - 2 * radius)) + radius;

person.posX = X;

person.posY = Y;

overlap = false;

for (j = 0; j < i; ++j) {

var person2 = people[j];

if (person2.Collision(person)) {

overlap = true;

break;

}

if (!overlap) break;

}

person.velX = Math.random() - 0.5;

person.velY = Math.random() - 0.5;

person.NormalizeVelocity();

person.velX *= speed;

person.velY *= speed;

if (i == 0) person.infected = true;

people.push(person);

}

Init();

function Advance() {

// for each from the population

for (i = 0; i < nrPeople; ++i) {

var person = people[i];

if (person.dead) continue;

// move

person.posX += person.velX * deltat;

person.posY += person.velY * deltat;

// collide / infect / cure

// first, walls collisions

if (person.posX <= radius) {

person.velX *= -1;

person.posX = radius;

}

else if (person.posX >= canvas.width - radius) {

person.velX *= -1;

person.posX = canvas.width - radius;

}

if (person.posY <= radius) {

person.velY *= -1;

person.posY = radius;

}

else if (person.posY >= canvas.height - radius) {

person.velY *= -1;

person.posY = canvas.height - radius;

}

// keep track of how long the infection lasts

if (person.infected)

person.infectedTime += deltat;

// collisions and infections between them

for (j = 0; j < i; ++j) {

var person2 = people[j];

if (person2.dead) continue;

if (person.Collision(person2)) {

person.Collide(person2);

if (person.infected && !person2.infected && !person2.cured) {

if (Math.random() < infectProb)

{

person2.infected = true;

++infections;

}

else if (person2.infected && !person.infected && !person.cured) {

if (Math.random() < infectProb)

{

person.infected = true;

++infections;

}

// cure

if (person.infected && person.infectedTime > cureTime) {

person.infected = false;

person.cured = true;

++cures;

}

// kill

if (person.infected && person.infectedTime > cureTime / 3.) {

if (Math.random() < deathProb * deltat / (cureTime * 2. / 3.)) {

person.dead = true;

++deaths;

}

// display

ctx.clearRect(0, 0, canvas.width, canvas.height);

for (i = 0; i < nrPeople; ++i) {

var person = people[i];

if (person.dead) continue;

ctx.beginPath();

ctx.arc(person.posX, person.posY, radius, 0, 2 * Math.PI, false);

ctx.fillStyle = person.infected ? 'red' : (person.cured ? 'green' : 'blue');

ctx.fill();

ctx.stroke();

}

canvasText.innerHTML = "Total: " + nrPeople + " Infected: " + infections + " Deaths: " + deaths + " Cured: " + cures + " Sick: " + (infections - cures - deaths);

ctxChart.beginPath();

ctxChart.moveTo(chartPosX * chartXScale, chart.height - chartValY * chartYScale);

++chartPosX;

chartValY = infections;

ctxChart.lineTo(chartPosX * chartXScale, chart.height - chartValY * chartYScale);

ctxChart.stroke();

if (infections - deaths == cures) // there is no more left to cure

{

ctxChart.clearRect(0, 0, chart.width, chart.height);

Init();

}

setInterval(Advance, 10);

})();

Conclusion

As usual, if you find mistakes or have suggestions, please point them out.

I’ll end up with a quote from Will the pandemic kill you? by Nobel Laureate Peter Doherty:

“What influenza also tells us is that, while viruses that spread via the respiratory route are the most likely cause of some future pandemic, only the most draconian and immediate restrictions on human travel are likely to limit the spread of infection, and then only briefly.”
“Such limitations are likely to be applied quickly if we are faced with a situation in which, for example, more than 30% of those affected develop severe or even fatal consequences. The more dangerous situation may be when mortality rates are in the 1-3% range, causing (ultimately) 70 million to 210 million deaths globally. Such an infection could “get away” before we realised what was happening.”

CodeProject

Epidemic

Introduction

Theory

Simulation

0

The code

Conclusion

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply

Introduction

Theory

Simulation

0

The code

Conclusion

Share this:

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply