После того, как заголовки о новом вирусе стали все больше заполнять ленты новостей, youtube и соц.сетей, стало интересно смоделировать ситуацию самостоятельно и посмотреть, на что влияют различные параметры. Сейчас уже много статей на Хабре на тему вируса, и даже появились сомнения, стоит ли выкладывать еще одну, но наверно хуже не будет. Многие вещи интуитивно и так понятны, но графически оно все же нагляднее, плюс можно менять параметры и сразу видеть результаты.
В конце текста приведен исходный код на языке Python, так что желающие могут поэкспериментировать самостоятельно.
Для тех, кому интересно, продолжение под катом.
В программной модели используется популяция размером 2000 человек. Число может быть любым, но слишком малые величины дают неточные результаты, слишком большие долго считаются. В симуляции каждый человек имеет свой «распорядок дня», массив на 24 значения с дискретом 1 час, включающий такие значения как «транспорт», «работа» или «прочее». Каждый человек может иметь двух возможных «родственников»: супруг и родитель. Еще два важных параметра, это возраст и статус заражения (здоров, болен, выздоровел и пр).
С точки зрения заражения, используется простая схема — первые 4-7 дней инфицированный человек может заражать других, затем он «заболевает» окончательно и попадает на карантин (никого не заражает), еще через несколько дней больной либо выздоравливает, либо нет (вероятность зависит от возраста), заражение второй раз невозможно. Минимальным шагом симуляции является 1 час. С каждым новым часом обновляется список людей в транспорте и на работе, пересчитывается количество инфицированных и выздоровевших, итоговые результаты выводятся в виде графика.
Посмотрим, что из всего этого можно получить.
Транспорт
В модели общественного транспорта используется 2 способа заражения:
— воздушно-капельный: если пассажир инфицирован, то он может заразить одного впереди, и одного сзади сидящего,
— контактный: если в салоне есть инфицированный пассажир, то заразиться может любой находящийся в автобусе, коснувшись например, ручки двери.
Таким образом, для каждого пассажира мы имеем сумму двух вероятностей, каждую из которых можно менять.
Разумеется, точных значений вероятности мы не знаем, поэтому посмотрим результаты для разных значений вероятности заражения: 4%, 1% и 0.25%. По вертикали графика процент числа зараженных, по горизонтали время в днях.
Виды графиков разумеется, условны, т.к. числа 4%, 1% и 0.25% взяты «из воздуха», но общий тренд видеть можно.
Работа
Тут все в принципе, то же самое. За тем исключением, что на работе человек проводит время, гораздо большее, чем в транспорте. Таким образом, уменьшение вероятности заражения приводит к более заметному падению значений на графике:
Здесь мы заодно можем наблюдать интересный парадокс — окончательное число переболевших мало зависит от вероятности заражения. График является более пологим, но одновременно, и более растянутым во времени. Это кажется странным, но если вдуматься, то это легко объяснимо: если человек ежедневно ездит в общественном транспорте, то уменьшение вероятности заражения вдвое приведет лишь к тому, что он заразится к примеру, не через 20, а через 40 поездок, но рано или поздно это все равно произойдет. С другой стороны, есть эффект, который врачи называют «сглаживанием кривой» (flattening the curve). К примеру, если в поликлинику или больницу одновременно обратится 20 человек, то им скорее всего окажут помощь и вылечат, если обратится 200 — то сделать это будет сложнее в силу недостатка лекарств, врачей и пр, если 2000 — надеюсь, читатели поняли идею…
На графике как раз хорошо видно, что в первом случае мы имеем резкий пик заболевших, с которым система здравоохранения может не справиться, в двух остальных случаях пики заметно меньше. Ну и в любом случае, в интересах человека уменьшить вероятность заражения по максимуму — время работает против болезни, разрабатываются новые методы лечения, появится вакцина и пр.
Семья
Здесь в модели используются два параметра вероятности. Первый это заражение между близкими родственниками (супруги, дети). Второй это заражение пожилых людей (родители). Если первый параметр достаточно сложно варьировать (в реальной жизни конечно, не в модели), т.к. сложно представить, мужа, жену и ребенка, живущих по отдельности, то посещение родителей минимизировать вполне можно, перейдя на общение по Скайпу, Вайберу и пр. И сделать это целесообразно, т.к. смертность у пожилых людей наиболее высока.
Приведу общий график где меняются оба параметра:
Общий принцип тут такой же, и подробно можно не комментировать.
Итоговый график
Выше мы меняли все параметры по отдельности. Наконец, настала пора объединить все параметры вместе и посмотреть результаты:
Здесь мы видим интересный результат. Только уменьшив все три параметра вероятности одновременно, мы можем выйти на тот режим в обществе, когда «коэффициент инфицирования» становится < 1, и эпидемия затухает без того, чтобы переболели все. Больные успевают выздороветь, не заразив новых больных. Но добиться этого достаточно сложно, нужно чтобы на всех этапах шанс распространиться у вируса был минимален.
Выводы
С медицинской точки зрения, выводов не будет — думаю, специализированных видео и статей уже достаточно.
С точки зрения математической, из графиков можно сделать несколько выводов:
— Есть вероятность что все-таки переболеют многие. С этим нужно просто смириться и принять как должное. Все же пневмония это не рак и не спид, по статистике (уже не модельной а реальной) для 80% даже не потребуется госпитализация.
— Несмотря на предыдущий пункт, минимизировать каждую из вероятностей заражения — можно и нужно, время работает на нас. ИТ-шники в этом плане имеют преимущество, т.к. им работать из дома вполне реально.
Как-то так. Будьте здоровы.
Для тех, кто захочет поэкспериментировать самостоятельно, исходник под спойлером. Все это делалось в свободное время, так что на красоту кода не претендую, ну и не гарантирую что нигде не ошибся.
# Coronavirus spreading simulation. For habr: https://habr.com/ru/post/492344/
import math
import random
import sys
import matplotlib.pyplot as plt
import numpy as np
from typing import List
population_count = 2000
# Person activity types
ACTIVITY_NONE = 0
ACTIVITY_TRANSPORT = 1
ACTIVITY_WORK = 2
ACTIVITY_SCHOOL = 3
# Capacity
BUS_CAPACITY = 40
NUM_BUSES = 24
WORKSPACE_CAPACITY = 40
SCHOOL_CAPACITY = 10
# Percentage of the population
BUS_USAGE_PERCENT = 50
# Population distribution
KIDS_PERCENT = 30
ADULTS_PERCENT = 40
RENTNERS_PERCENT = 30
WORKERS_PERCENT = 60
# Person status
NORMAL = 1
INFECTED = 2
SICK = 3
IMMUNE = 4
DEAD = 5
class Person:
def __init__(self):
# Family
self.parent = None
self.spouse = None
# Day activities
self.activities = 24*[ACTIVITY_NONE]
self.activity = ACTIVITY_NONE
self.age = 0
# Infection status
self.status = NORMAL
self.day_infected = -1
def make_kid(self):
t_start = random.randint(7,9)
self.activities[t_start+1:t_start+9+1] = 9*[ACTIVITY_SCHOOL]
self.age = random.randint(3,17)
return self
def make_adult(self):
t_start = random.randint(7,9)
use_bus = BUS_USAGE_PERCENT > random.randint(0, 100)
self.activities[t_start] = ACTIVITY_TRANSPORT if use_bus else ACTIVITY_NONE
self.activities[t_start+1:t_start+9+1] = 9*[ACTIVITY_WORK]
self.activities[t_start+9+1] = ACTIVITY_TRANSPORT if use_bus else ACTIVITY_NONE
self.age = random.randint(18, 68)
return self
def make_rentner(self):
self.activities[random.randint(11, 18)] = ACTIVITY_TRANSPORT
self.activities[random.randint(11, 18)] = ACTIVITY_TRANSPORT
self.age = random.randint(68, 88)
return self
def is_kid(self):
return self.age >= 3 and self.age <= 17
def is_adult(self):
return self.age >= 18 and self.age <= 68
def is_infected(self):
return self.status == INFECTED
def is_athome(self, day_hour):
return self.activities[day_hour] != ACTIVITY_WORK and self.activities[day_hour] != ACTIVITY_WORK
def death_probability(self):
# https://www.worldometers.info/coronavirus/coronavirus-age-sex-demographics/
if self.age < 10:
return 0.0
if self.age < 20:
return 0.2
if self.age < 40:
return 0.4
if self.age < 50:
return 0.8
if self.age < 75:
return 5.0
if self.age < 80:
return 8.0
return 14.0
def infect(self, probability: float, day_num: int):
if random.random()*100.0 < probability and self.status == NORMAL:
self.status = INFECTED
self.day_infected = day_num
def update_infection(self, hour: int):
day_num = hour // 24
n_days = day_num - self.day_infected
if self.status == INFECTED and n_days >= random.randint(4, 8):
self.status = SICK
elif self.status == SICK and n_days >= random.randint(10, 20):
if random.random() < 0.01 * self.death_probability():
self.status = DEAD
else:
self.status = IMMUNE
def print(self):
print(f"Person: {self.activities}, Age: {self.age}, Status: {self.status}")
class Bus:
def __init__(self):
self.capacity = BUS_CAPACITY
self.passengers:List[Person] = BUS_CAPACITY*[None]
def add_passenger(self, passenger: Person):
# Get indexes of free seats
places = []
for p in range(self.capacity):
if self.passengers[p] is None:
places.append(p)
# Passenger use a random seat
if len(places) > 0:
new_index = random.randint(0, len(places) - 1)
# print(" ", places, places[new_index])
self.passengers[places[new_index]] = passenger
def del_passengers(self, day_hour:int):
# Passenges leaves the bus if their travel time is over
for p in range(self.capacity):
self.passengers[p] = None
def free_seats(self):
places = 0
for p in range(self.capacity):
if self.passengers[p] is None:
places += 1
return places
def update_infection(self, hour: int, probability: float):
day_num = hour//24
is_somebody_infected = False
for p in range(len(self.passengers)):
if self.passengers[p] is not None and self.passengers[p].is_infected():
is_somebody_infected = True
# Spread to the nearest seats (air)
if p > 0 and self.passengers[p-1] is not None:
self.passengers[p - 1].infect(probability=probability, day_num=day_num)
if p < len(self.passengers) - 1 and self.passengers[p+1] is not None:
self.passengers[p + 1].infect(probability=probability, day_num=day_num)
# Spread to all passengers (hands touch)
if is_somebody_infected:
for person in self.passengers:
if person is not None:
person.infect(probability=probability, day_num=day_num)
def print(self):
s = ''
for p in self.passengers:
s += f"{p.status if p else '-'} "
print(f"Bus: {self.free_seats()}/{len(self.passengers)} seats, Passengers: {s}")
class Work:
def __init__(self):
self.capacity = WORKSPACE_CAPACITY
self.people = []
def free_places(self):
return self.capacity - len(self.people)
def add_person(self, person: Person):
self.people.append(person)
def update_infection(self, hour: int, probability: float):
day_num = hour // 24
day_hour = hour % 24
# Spread to the nearest places (air)
is_somebody_infected = False
for p, person in enumerate(self.people):
if person.activities[day_hour] == ACTIVITY_WORK and person.is_infected():
is_somebody_infected = True
if p > 0:
self.people[p - 1].infect(probability=probability, day_num=day_num)
if p < len(self.people) - 1:
self.people[p + 1].infect(probability=probability, day_num=day_num)
# Spread to all people (hands touch)
if is_somebody_infected:
for person in self.people:
person.infect(probability=probability, day_num=day_num)
def print(self):
s = ''
for p in self.people:
s += f"{p.status} "
print(f"Workplace: {len(self.people)} people: {s}")
class Population:
def __init__(self):
self.people: List[Person] = []
self.transport: List[Bus] = []
self.workspaces: List[Work] = []
self.probability_infection_transport = 2.0
self.probability_infection_work = 2.0
self.probability_infection_family = 2.0
self.probability_infection_grandfamily = 2.0
def create_population(self):
n_kids = KIDS_PERCENT * population_count // 100
n_adults = ADULTS_PERCENT * population_count // 100
n_rentners = RENTNERS_PERCENT * population_count // 100
print(f"Creating the population: {population_count} people, kids: {n_kids}, adults: {n_adults}, rentners: {n_rentners}")
# Add population
self.people.clear()
adults, rentners, kids = [], [], []
for _ in range(n_adults):
person = Person().make_adult()
self.people.append(person)
adults.append(person)
# Add person to a workspace
for work in self.workspaces:
if work.free_places() > 0:
work.add_person(person)
break
for _ in range(n_kids):
kid = Person().make_kid()
kids.append(kid)
self.people.append(kid)
for _ in range(n_rentners):
rentner = Person().make_rentner()
self.people.append(rentner)
rentners.append(rentner)
# Add parents
for kid in kids:
kid.parent = random.choice(adults)
for person in adults:
if random.random() > 0.5:
person.parent = random.choice(rentners)
# Add spouses: half of people are married
for person in adults:
if random.random() > 0.5:
person.spouse = random.choice(adults)
for person in rentners:
if random.random() > 0.75:
person.spouse = random.choice(rentners)
# Add one infected random person
random.choice(adults).status = INFECTED
def create_transport(self):
print(f"Creating the transport: {population_count*BUS_USAGE_PERCENT//100} people using transport, {NUM_BUSES} buses total")
for _ in range(NUM_BUSES):
self.transport.append(Bus())
def create_workspaces(self):
n_adults = ADULTS_PERCENT * population_count // 100
n_workspaces = math.ceil(n_adults / (WORKSPACE_CAPACITY))
print(f"Creating {n_workspaces} workspace, each of {WORKSPACE_CAPACITY} people")
for _ in range(n_workspaces):
self.workspaces.append(Work())
def update_transport(self, hour):
day_hour = hour % 24
# Passengers are leaving the transport
for bus in self.transport:
bus.del_passengers(day_hour)
# Passengers are entering the transport
for person in self.people:
if person.activities[day_hour] == ACTIVITY_TRANSPORT:
# Find random bus for the passenger
# buses = []
# for bus in self.transport:
# if bus.free_seats() > 0:
# buses.append(bus)
buses = list(filter(lambda bus: bus.free_seats() > 0, self.transport))
if len(buses) > 0:
bus = random.choice(buses)
bus.add_passenger(person)
def update_infection(self, hour):
day_num = hour // 24
day_hour = hour % 24
# Transport
for bus in self.transport:
bus.update_infection(hour, probability=self.probability_infection_transport)
# Workspace
for work in self.workspaces:
work.update_infection(hour, probability=self.probability_infection_work)
# Family
for person in self.people:
# Person <=> Parent
if person.is_infected() and person.is_athome(day_hour) and person.parent is not None and person.parent.is_athome(day_hour):
prob = self.probability_infection_family if person.parent.age < 68 else self.probability_infection_grandfamily
person.parent.infect(probability=prob, day_num=day_num)
elif person.parent is not None and person.is_athome(day_hour) and person.parent.is_infected() and person.parent.is_athome(day_hour):
person.infect(probability=self.probability_infection_family, day_num=day_num)
# Person <=> Spouse
if person.is_infected() and person.is_athome(day_hour) and person.spouse is not None and person.spouse.is_athome(day_hour):
person.spouse.infect(probability=self.probability_infection_family, day_num=day_num)
elif person.spouse is not None and person.spouse.is_athome(day_hour) and person.spouse.is_infected() and person.spouse.is_athome(day_hour):
person.infect(probability=self.probability_infection_family, day_num=day_num)
# Update health status for infected people
for person in self.people:
person.update_infection(hour)
def get_status(self, hour):
normal, infected, sick, dead, immune = 0, 0, 0, 0, 0
for person in self.people:
if person.status == INFECTED:
infected += 1
elif person.status == IMMUNE:
immune += 1
elif person.status == SICK:
sick += 1
elif person.status == DEAD:
dead += 1
elif person.status == NORMAL:
normal += 1
print(f"Day {hour//24}: Healthy: {normal}, Infected: {infected}, Sick: {sick}, Dead: {dead}, Immune: {immune}")
return normal, infected, sick, dead, immune
def test_buses(self):
p = Person()
p.make_adult()
b = Bus()
b.print()
b.add_passenger(p)
b.print()
b.add_passenger(Person())
b.print()
b.add_passenger(Person())
b.print()
if __name__ == "__main__":
print("Application started")
country1 = Population()
country1.probability_infection_transport = 2.0
country1.probability_infection_work = 2.0
country1.probability_infection_family = 2.0
country1.probability_infection_grandfamily = 2.0
country2 = Population()
country2.probability_infection_transport = 1.0
country2.probability_infection_work = 1.0
country2.probability_infection_family = 1.0
country2.probability_infection_grandfamily = 1.0
country3 = Population()
country3.probability_infection_transport = 0.25
country3.probability_infection_work = 0.25
country3.probability_infection_family = 0.25
country3.probability_infection_grandfamily = 0.0
country1.create_transport()
country1.create_workspaces()
country1.create_population()
country2.create_transport()
country2.create_workspaces()
country2.create_population()
country3.create_transport()
country3.create_workspaces()
country3.create_population()
print()
print("Starting the simulationn")
days_list = []
normal_list1, infected_list1, immune_list1, dead_list1 = [], [], [], []
normal_list2, infected_list2, immune_list2, dead_list2 = [], [], [], []
normal_list3, infected_list3, immune_list3, dead_list3 = [], [], [], []
hour = 0
for day in range(1, 96):
for hr in range(24):
country1.update_transport(hour)
country1.update_infection(hour)
country2.update_transport(hour)
country2.update_infection(hour)
country3.update_transport(hour)
country3.update_infection(hour)
# if hr >= 6 and hr <= 20:
# for bus in country1.transport:
# bus.print()
# for work in country1.workspaces:
# work.print()
# print()
hour += 1
days_list.append(day)
normal, infected, _, dead, immune = country1.get_status(hour)
normal_list1.append(normal*100/population_count)
infected_list1.append(infected*100/population_count)
immune_list1.append(immune*100/population_count)
dead_list1.append(dead*100/population_count)
normal, infected, _, dead, immune = country2.get_status(hour)
normal_list2.append(normal*100/population_count)
infected_list2.append(infected*100/population_count)
immune_list2.append(immune*100/population_count)
dead_list2.append(dead*100/population_count)
normal, infected, _, dead, immune = country3.get_status(hour)
normal_list3.append(normal*100/population_count)
infected_list3.append(infected*100/population_count)
immune_list3.append(immune*100/population_count)
dead_list3.append(dead*100/population_count)
plt.rcParams["figure.figsize"] = (12, 4)
plt.plot(days_list, infected_list1, label="Infected, prob. 2%")
plt.plot(days_list, infected_list2, label="Infected prob. 1%")
plt.plot(days_list, infected_list3, label="Infected prob. 0.25%")
# plt.plot(days_list, dead_list1, label="Dead 2%")
# plt.plot(days_list, dead_list2, label="Dead 1%")
# plt.plot(days_list, dead_list3, label="Dead 0.25%")
plt.ylim(-5, 50)
plt.xticks(np.arange(min(days_list), max(days_list) + 1, 30.0))
plt.legend()
plt.show()
print("Done")
print()
