# 7. Differential Equations

This step illustrates how to use differential equations.

## Formulation

We are interested in the spreading of the disease inside the buildings. In order to model it, we will use differential equations. So, we will need to:

- Add two global variables to define the building epidemic properties (
) and numerical integration parameter (`beta`

).`h`

- Add new variables for the buildings (
,`I`

,`S`

,`T`

,`t`

) to manage epidemic;`I_to1`

- Define differential equations for disease spreading inside buildings.
- Add one behavior for buildings for the spreading of the disease.

## Model Definition

### global variables

We define two new global variables used in the disease spreading dynamic inside the buildings: (i) `beta`

is the contamination rate, and `h`

is the integration step (used in the `solve`

statement).

```
global {
...
float beta <- 0.01;
float h <- 0.1;
...
}
```

### building

In order to define the disease spread dynamics, we define several variables that will be used by the differential equations:

: float, number of people infected in the building.`I`

: float, number of people not infected in the building.`S`

: float, the total number of people in the building.`T`

: float, the current time of the equation system integration.`t`

: float, the remaining number of people infected (float number lower between 0 and 1 according to the differential equations).`I_to1`

```
species building {
...
float I;
float S;
float T;
float t;
float I_to1;
...
}
```

Then, we define the differential equations system that will be used for the disease spreading dynamic. Note that to define a differential equation system we use the block ** equation** + name. These equations are the classic ones used by SI mathematical models.

```
species building {
....
equation SI{
diff(S,t) = (- beta * S * I / T) ;
diff(I,t) = ( beta * S * I / T) ;
}
...
}
```

At last, we define a new reflex for the building called ** epidemic** that will be activated only when there is someone inside the building. This reflex first computes the number of people inside the building (

`T`

), then the number of not infected people (`S`

) and finally the number of infected ones (`I`

).If there is at least one people infected and one people not infected, the differential equations is integrated (according to the integration step value `h`

) with the method Runge-Kutta 4 to compute the new value of infected people. We then sum the old value of `I_to1`

with the number of people newly infected (this value is a float and not an integer). Finally, we cast this value as an integer, ask the corresponding number of not infected people to become infected, and decrement this integer value to `I\_to1`

.

```
species building {
...
reflex epidemic when: nb_total > 0 {
T <- float(nb_total);
S <- float(nb_total - nb_infected);
I <- T - S;
float I0 <- I;
if (I > 0 and S > 0) {
solve SI method: #rk4 step_size: h;
I_to1 <- I_to1 + (I - I0);
int I_int <- min([int(S), int(I_to1)]);
I_to1 <- I_to1 - I_int;
ask (I_int among (people_in_building where (!each.is_infected))) {
is_infected <- true;
}
}
}
...
}
```

## Complete Model

```
/**
* Name: Differential Equation
* Author: GAMA Team
* Description: 7th part of the tutorial : Incremental Model
* Tags: tutorial, chart, graph, 3d, light, multi-Level, equation
*/
model model7
global {
int nb_people <- 500;
float agent_speed <- 5.0 #km / #h;
float step <- 1 #minutes;
float infection_distance <- 2.0 #m;
float proba_infection <- 0.05;
int nb_infected_init <- 5;
file roads_shapefile <- file("../includes/road.shp");
file buildings_shapefile <- file("../includes/building.shp");
geometry shape <- envelope(roads_shapefile);
graph road_network;
float staying_coeff update: 10.0 ^ (1 + min([abs(current_date.hour - 9), abs(current_date.hour - 12), abs(current_date.hour - 18)]));
float beta <- 0.01;
float h <- 0.1;
list<people_in_building> list_people_in_buildings update: (building accumulate each.people_in_building);
int nb_people_infected <- nb_infected_init update: (people + list_people_in_buildings) count (each.is_infected);
int nb_people_not_infected <- nb_people - nb_infected_init update: nb_people - nb_people_infected;
bool is_night <- true update: current_date.hour < 7 or current_date.hour > 20;
float infected_rate update: nb_people_infected / nb_people;
init {
create road from: roads_shapefile;
road_network <- as_edge_graph(road);
create building from: buildings_shapefile;
create people number: nb_people {
speed <- agent_speed;
building bd <- one_of(building);
location <- any_location_in(bd);
}
ask nb_infected_init among people {
is_infected <- true;
}
}
reflex end_simulation when: infected_rate = 1.0 {
do pause;
}
}
species people skills: [moving] {
bool is_infected <- false;
point target;
int staying_counter;
reflex move when: target != nil {
do goto target: target on: road_network;
if (location = target) {
target <- any_location_in(one_of(building));
target <- nil;
staying_counter <- 0;
}
}
reflex infect when: is_infected {
ask people at_distance infection_distance {
if flip(proba_infection) {
is_infected <- true;
}
}
}
aspect circle {
draw circle(5) color: is_infected ? #red : #green;
}
aspect sphere3D {
draw sphere(3) at: {location.x, location.y, location.z + 3} color: is_infected ? #red : #green;
}
}
species road {
geometry display_shape <- shape + 2.0;
aspect default {
draw display_shape color: #black depth: 3.0;
}
}
species building {
int nb_infected <- 0 update: self.people_in_building count each.is_infected;
int nb_total <- 0 update: length(self.people_in_building);
float height <- rnd(10 #m, 20 #m);
float I;
float S;
float T;
float t;
float I_to1;
species people_in_building parent: people schedules: [] { }
reflex let_people_leave {
ask self.people_in_building {
staying_counter <- staying_counter + 1;
}
release people_in_building where (flip((each).staying_counter / staying_coeff)) as: people in: world {
target <- any_location_in(one_of(building));
}
}
reflex let_people_enter {
capture (people inside self) where (each.target = nil) as: people_in_building;
}
equation SI {
diff(S, t) = (-beta * S * I / T);
diff(I, t) = (beta * S * I / T);
}
reflex epidemic when: nb_total > 0 {
T <- float(nb_total);
S <- float(nb_total - nb_infected);
I <- T - S;
float I0 <- I;
if (I > 0 and S > 0) {
solve SI method: #rk4 step_size: h;
I_to1 <- I_to1 + (I - I0);
int I_int <- min([int(S), int(I_to1)]);
I_to1 <- I_to1 - I_int;
ask (I_int among (people_in_building where (!each.is_infected))) {
is_infected <- true;
}
}
}
aspect default {
draw shape color: nb_total = 0 ? #gray : (float(nb_infected) / nb_total > 0.5 ? #red : #green) border: #black depth: height;
}
}
experiment main_experiment type: gui {
parameter "Infection distance" var: infection_distance;
parameter "Proba infection" var: proba_infection min: 0.0 max: 1.0;
parameter "Nb people infected at init" var: nb_infected_init;
output {
monitor "Current hour" value: current_date.hour;
monitor "Infected people rate" value: infected_rate;
display map_3D type: opengl {
light 1 color: (is_night ? 50 : 255);
image "../includes/soil.jpg";
species road;
species people aspect: sphere3D;
species building transparency: 0.5;
}
display chart refresh: every(10 #cycles) {
chart "Disease spreading" type: series {
data "susceptible" value: nb_people_not_infected color: #green marker: false;
data "infected" value: nb_people_infected color: #red marker: false;
}
}
}
}
```