Симуляция распространения коронавируса

Привет Хабр.

После того, как заголовки о новом вирусе стали все больше заполнять ленты новостей, 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% даже не потребуется госпитализация.
— Несмотря на предыдущий пункт, минимизировать каждую из вероятностей заражения — можно и нужно, время работает на нас. ИТ-шники в этом плане имеют преимущество, т.к. им работать из дома вполне реально.

Как-то так. Будьте здоровы.

Для тех, кто захочет поэкспериментировать самостоятельно, исходник под спойлером. Все это делалось в свободное время, так что на красоту кода не претендую, ну и не гарантирую что нигде не ошибся.

covid.py

# 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()

 

Источник

, , , , , ,

Читайте также

Меню