domain: info sdv

text:
Ayan Chaudhury
email: ayan.chaudhury@inria.fr

Frédéric Boudon
email: frederic.boudon@cirad.fr

Christophe Godin
email: christophe.godin@inria.fr

3D Plant Phenotyping: All You Need is Labelled Point Cloud Data

Keywords: Procedural Model, L-System, Virtual Plant, Point Cloud, Labeled Synthetic Data

In the realm of modern digital phenotyping technological advancements, the demand of annotated datasets is increasing for either training machine learning algorithms or evaluating 3D phenotyping systems. While a few 2D datasets have been proposed in the community in last few years, very little attention has been paid to the construction of annotated 3D point cloud datasets. There are several challenges associated with the creation of such annotated datasets. Acquiring the data requires instruments having good precision and accuracy levels. Reconstruction of full 3D model from multiple views is a challenging task considering plant architecture complexity and plasticity, as well as occlusion and missing data problems. In addition, manual annotation of the data is a cumbersome task that cannot easily be automated. In this context, the design of synthetic datasets can play an important role. In this paper, we propose an idea of automatic generation of synthetic point cloud data using virtual plant models. Our approach leverages the strength of the classical procedural approach (like L-systems) to generate the virtual models of plants, and then perform point sampling on the surface of the models. By applying stochasticity in the procedural model, we are able to generate large number of diverse plant models and the corresponding point cloud data in a fully automatic manner. The goal of this paper is to present a general strategy to generate annotated 3D point cloud datasets from virtual models. The code (along with some generated point cloud models) are available at: https://gitlab.inria.fr/mosaic/publications/lpy2pc.

Introduction

With the recent breakthrough of the advancements on plant phenotyping research, botanical and agronomic studies have been revolutionized to enter in a new stage of development [START_REF] Tardieu | Plant phenomics, from sensors to knowledge[END_REF]. In recent years, 3D computer vision based systems have been widely developed in different types of phenotyping applications including robotic branch pruning [START_REF] Chattopadhyay | Measuring and modeling apple trees using time-of-flight data for automation of dormant pruning applications[END_REF], automated growth analysis [START_REF] Chaudhury | Computer vision based autonomous robotic system for 3d plant growth measurement[END_REF][START_REF] Chaudhury | Machine vision system for 3d plant phenotyping[END_REF], agricultural automation tasks [START_REF] Vázquez-Arellano | 3-d imaging systems for agricultural applications -a review[END_REF], etc. Point cloud based approach has been an integral part of machine learning aspects of plant phenotyping [START_REF] Ziamtsov | Machine learning approaches to improve three basic plant phenotyping tasks using three-dimensional point clouds[END_REF]. However, many of these techniques require ground truth data for training the algorithms and for validation of the results in a quantitative manner. For example, methods performing plant organ segmentation [START_REF] Paulus | Surface feature based classification of plant organs from 3d laserscanned point clouds for plant phenotyping[END_REF][START_REF] Li | Analyzing growing plants from 4d point cloud data[END_REF] require labeled point cloud data for quantifying the segmentation accuracy. Similarly, in-field 3D segmentation [START_REF] Sodhi | In-field segmentation and identification of plant structures using 3d imaging[END_REF] or plant structure classification tasks [START_REF] Dey | Classification of plant structures from uncalibrated image sequences[END_REF] also demand annotated datasets for validation. In addition, to exploit the recent success of deep learning technologies for point cloud data [START_REF] Qi | Pointnet: Deep learning on point sets for 3d classification and segmentation[END_REF] in agricultural applications, large amount of annotated data are in extreme demand. Unfortunately at the current moment, there is a scarcity of available annotated point cloud dataset in the plant phenotyping community. Only a handful of datasets of certain plant species are proposed in the literature. Among the currently available 3D datasets, Cruz et al. [START_REF] Cruz | Multimodality imagery database for plant phenotyping[END_REF] created a multimodal database of top view depth images of Arabidopsis and bean plants. Wen et al. [44] proposed a database of 3D geometric plant models which are based on the measurements of in-situ morphological data. Different types of species in different growth stages are modelled in the framework. However the data is limited to the organ scale of the plants, and not on the full 3D model. Bernotas et al. [START_REF] Bernotas | A photometric stereo-based 3d imaging system using computer vision and deep learning for tracking plant growth[END_REF] created a photometric stereo dataset along with manually annotated leaf masks of Arabidopsis rosettes. The Komatsuna plant dataset [START_REF] Uchiyama | An easy-to-setup 3d phenotyping platform for komatsuna dataset[END_REF] contains RGB and depth images with annotated leaf labels. However, none of the above mentioned datasets are built on full 3D point cloud model of plants. Recently, the ROSE-X dataset [START_REF] Dutagaci | Rose-x: an annotated data set for evaluation of 3d plant organ segmentation methods[END_REF] is created with 3D models of rosebush plants acquired through X-ray tomography method. Manual annotation is performed to label voxels and point clouds belonging to different organs of the plant. For remote sensing applications of vegetation studies, Wang [START_REF] Wang | Unsupervised semantic and instance segmentation of forest point clouds[END_REF] simulated a synthetic Terrestrial Laser Scanning (TLS) data of forest scenes. The generated point cloud contains simulated data of large trees having structurally complex crowns with labels belonging to leaf and wood class.

In contrast to the scarcity of 3D point cloud datasets, there has been quite a few 2D datasets available for public use, as well as software tool to access different types of existing solutions according to the experimental need [START_REF] Lobet | An online database for plant image analysis software tools[END_REF]. These have revealed extremely useful for training machine learning algorithms and performance validations tasks for 2D image based plant phenotyping applications. The dataset proposed in [START_REF] Minervini | Finely-grained annotated datasets for image-based plant phenotyping[END_REF] is one of the first kind, where annotated top-view color images of rosette plants at different growth stages are available. Shadrin et al. [START_REF] Shadrin | Instance segmentation for assessment of plant growth dynamics in artificial soilless conditions[END_REF] proposed a manually annotated dataset of raw salad images acquired over a period of several days for plant growth dynamics assessment. Among other recent datasets, the Oil Radish dataset [START_REF] Mortensen | The oil radish growth dataset for semantic segmentation and yield estimation[END_REF] and Grass clover dataset [START_REF] Skovsen | The grassclover image dataset for semantic and hierarchical species understanding in agriculture[END_REF] have been proposed. In last few years, leaf segmentation and counting of rosette plants have been an active area of interest. Deep learning techniques have been very successful in solving this type of problem, which demand lot of annotated data for training the network. In order to overcome the scarcity of available datasets, different types of data augmentation strategies have been proposed to generate realistic synthetic data. Giuffrida et al. [START_REF] Giuffrida | Arigan: Synthetic arabidopsis plants using generative adversarial network[END_REF] proposed a Generative Adversarial Network (GAN) model to artificially create Arabidopsis rosette images. The model allows the user to create plants with desired number of leaves in a realistic manner. With a similar type of motivation, Ward et al. [START_REF] Ward | Deep leaf segmentation using synthetic data[END_REF] generated synthetic leaf images inspired by domain randomization strategy using a Mask-RCNN deep learning model. Kuznichov et al. [START_REF] Kuznichov | Data augmentation for leaf segmentation and counting tasks in rosette plants[END_REF] proposed a data augmentation strategy which preserves the geometric structure of the generated leaves with the real data. Different types of synthetic images are generated by applying some heuristic rules on the leaf shape and growth and recently, parametric L-systems have been used to perform data augmentation in generating synthetic leaf models [START_REF] Ubbens | The use of plant models in deep learning: an application to leaf counting in rosette plants[END_REF]. Various leaf databases have also been proposed [START_REF] Soderkvist | Computer vision classification of leaves from swedish trees[END_REF][START_REF] Wu | A leaf recognition algorithm for plant classification using probabilistic neural network[END_REF][START_REF] Kumar | Leafsnap: A computer vision system for automatic plant species identification[END_REF] which are useful for species identification and leaf architecture analysis [START_REF] Wilf | Computer vision cracks the leaf code[END_REF][START_REF] Wang | Can walking and measuring along chord bunches better describe leaf shapes?[END_REF][START_REF] Chaudhury | Plant species identification from occluded leaf images[END_REF]. Some of the difficulties associated with creating labeled 3D point cloud datasets include the expensive scanning devices, cumbersome and error prone manual annotation of large number of points, and reconstructing full 3D structure of complex plants in the presence of occlusion and noise. Also, growing varieties of plant species and acquiring data at different growth stages require lot of manpower and constant monitoring of the plants. In this regard, synthetic (or virtual) plant models can play an important role. Recently Ward et al. [START_REF] Ward | Scalable learning for bridging the species gap in imagebased plant phenotyping[END_REF] demonstrated that synthetic data can be extremely useful in training deep learning models, which can be reliably used for measurements of different types of phenotypic traits in the general case. Modelling the geometry of complex plant structures have been a center of attention for mathematicians and biologists for decades [START_REF] Lindenmayer | The algorithmic beauty of plants[END_REF][START_REF] Godin | Functional-structural plant modelling[END_REF][START_REF] Godin | Plant architecture modelling -virtual plants and complex systems[END_REF][START_REF] Prusinkiewicz | Computational models of plant development and form[END_REF]. Virtual models can be extremely useful in agricultural studies [START_REF] Evers | Modeling branching in cereals[END_REF][START_REF] Buck-Sorlin | Meeting present and future challenges in sustainable horticulture using virtual plants[END_REF], as well as have great potential for mechanical simulation of plant behavior [START_REF] Bucksch | Morphological plant modeling: unleashing geometric and topological potential within the plant sciences[END_REF]. In generating the virtual plant models, L-system based modelling technique has been quite successful [START_REF] Prusinkiewicz | The use of positional information in the modeling of plants[END_REF]. Different types of platforms have been developed to simulate the L-system rules, e.g., L+C modelling language [START_REF] Karwowski | Design and implementation of the l+c modeling language[END_REF], L-Py framework [START_REF] Boudon | L-py: an lsystem simulation framework for modeling plant architecture development based on a dynamic language[END_REF], etc. Functional Structural Plant Models (FSPM) [START_REF] Godin | Functional-structural plant modelling[END_REF] have emerged as powerful tool to construct 3D models of plant functioning and growth. These models can play an important role in mechanical simulation in crop science research [START_REF] Vos | Functional-structural plant modelling: a new versatile tool in crop science[END_REF]. Although generation of synthetic plant models have been a well studied area, the virtual plant models have not been exploited in generating synthetic 3D point cloud data for plant phenotyping applications. Compared to the synthetic artery data simulation software like Vascusynth [START_REF] Hamarneh | Vascusynth: Simulating vascular trees for generating volumetric image data with ground-truth segmentation and tree analysis[END_REF] for medical imaging research, there is no such tool available in the plant phenotyping community.

Inspired by the necessity of synthetic plant point cloud datasets for 3D plant phenotyping applications, we propose a general approach to sample points on the surface of virtual plant models to generate annotated 3D point cloud data in a fully automatic manner. The framework is general in the sense that it can be plugged into any procedural plant model, and surface point cloud can be sampled from the virtual plant. Another motivation of the work comes from a computer graphics perspective. Among different types of representations of geometric shapes in computer graphics, two types are widely used: triangular mesh and point based representation [START_REF] Kobbelt | A survey of point-based techniques in computer graphics[END_REF]. These two types are complementary to each other in terms of information richness and model complexity, and our framework builds a bridge between the two paradigms for the varieties of applications of the geometric plant models in phenotyping problems.

Initially we represent the virtual model as a collection of geometric primitives, each having their own label belonging to different organs of a plant. Each primitive is then tessellated into triangles to obtain a basic coarse mesh model. By a random sampling strategy using the idea of barycentric coordinates, we are able to generate uniformly distributed point cloud on the surface of the virtual model. Points residing inside any primitive are removed by a geometric filtering technique. We also model some basic characteristics of real world scanning devices by incorporating different sampling strategies, which can be adopted according to the need of the application. We demonstrate the efficiency of the approach by generating realistic synthetic point cloud model of Arabidopsis plants, as well as plants with different types of architecture. These datasets are generated by the stochasticity inherent to the procedural model and by a second level of stochasticity at the level of point cloud sampling. We also show some examples on other synthetic plants having different types of geometry structures.

In the next section we briefly describe how a virtual plant model should be specified to allow our algorithm to construct annotated point clouds from it. Then the basic random point cloud sampling strategy is explained, followed by the insideness testing process. Then we discuss the sampling strategy to obtain device specific point cloud models, and experimental results are demonstrated. Finally the scope of the work and directions for further improvements are discussed.

Virtual Plant Models

Virtual plants are computer generated synthetic models mimicking the development and growth of real plants. The models are typically developed by the recursive rules of L-systems [START_REF] Lindenmayer | The algorithmic beauty of plants[END_REF]. One of the main attribute of L-system based plant models is that, complex plant structures can be generated by successive application of simple recursive rules. The recursive rules help in modelling natural phenomena like plant growth, creation and destruction of plant organs over time, change of shape (e.g. bending, spiralling), etc. This helps us to model a desired plant species in different growth stages by incorporating biological knowledge about the plant in the L-system rules.

In a L-system framework, a plant is defined by a string of symbols, called L-string, representing the different organs of the plant, A for Apex, I for internode, L for leaf for example. Each symbol is given a set of arguments describing variables attached to the corresponding organ. For example an internode can be modeled as a cylinder of height h and radius r, a leaf can be modeled as a polygon of size s, an apex can have an age t, etc. Square brackets are used to delineate branches within the L-string. An initial L-string, called the axiom, defines the initial state of the plant, before growth starts. For example, the plant may initially be composed of an apex A(t). This is represented as: Axiom: A(0) # meaning that the L-string initially contains a single apex at time 0

A typical L-system model, then consists of a set of rules that make it possible to evolve the initial L-string. For instance the following rule typically says that at each time step dt the apex produces an internode I, a branching apex A and a leaf L, and an apical apex A:

A(t) --> I(0.01,0.1) [+(angle)[L(0.1)]A(t+dt)] A(t+dt)
where 'angle' is a variable (that can be a constant for example). This rule can be complemented by rules that express how internodes and leaf should change in time:

I(r,h) --> I(r',h') with r' = r + a*dt and h' = r' + b*dt L(s) --> L(s') with s' = s + c*dt
In general, these so-called production rules do not describe how the symbols are associated with geometric models. For this, a way to interpret geometrically the symbols of the L-string should be defined. This is done using a new set of rules called interpretation rules. Here they could for instance look like, e.g.:

A(t) --> Sphere(radius(t))
with radius(t) = 0.1 (constant radius in time) I(r,h) --> Cylinder(r,h) L(s) --> Polygon(size)

To use our algorithm, one would need to attach a unique annotation for each class of organ that one would like to recognize in the application. This is done by adding an extra 'label' argument to the symbols of interest.

A(t) --> I(0.01,0.1,8) [+(angle)[L(0.1,9)]A(t+dt)] A(t+dt) I(r,h,label_I) --> I(r',h',Label_I) with r' = r + a*dt and h' = r' + b*dt L(size,label_L) --> L(s',label_L) with s' = s + c*dt
where label=8 is attached to internodes and label=9 to leaves. The interpretation rules should be modified accordingly:

A(t) --> Sphere(radius(t))
with radius(t) = 0.1 (constant radius in time) I(r,h,label_I) --> Cylinder(r,h,label_I) L(size,label_L) --> Polygon(size,label_L)

In this way, the geometric models will be tagged with a label corresponding to their category (here Internode or Leaf), and will make it possible to identify the triangles (which are used for rendering the model) that belong respectively to internodes and leaves in the 3D scene. Basically we transfer the label of an organ to the triangles while performing the tessellation, as discussed next.

Surface Point Cloud Generation

After creating the virtual model, we aim at generating points on the surface of the model where each point has a label of the plant organ it belongs to. In order to achieve this, we follow a series of steps. Given a parametric L-system model with unique label associated with each organ of the plant, we extract the set of primitives (along with their label of the organ they belong to in the procedural model) which are used to display the virtual plant model as a 3D object. In general, the following types of primitives are typically used in constructing the 3D object: cylinder, frustum, sphere, Beźier curve, NURBS patch, etc. Each primitive is then locally tessellated into a set of triangles (these triangles are used to render the 3D object), and the triangles are assigned the label of the primitive. Finally a global list of labelled triangle vertices are created to obtain a coarse mesh model of the 3D structure of the plant. Now the idea is to sample a dense set of labelled point cloud on the triangle surfaces, with a controlled density to create a realistic point cloud model of the plant. In the next section, we describe a simple and effective random point sampling strategy that can generate desired number of labelled points on the surface of the model.

Random Point Sampling

In order to generate labelled points on the surface of the virtual model, we propose two types of sampling strategies. The first type of sampling aims at generating uniform density of points on the model surface, and the second type of sampling is motivated by the effect of real world scanning devices, for which point density depends on various sensor-related parameters. In this section, we present the basic sampling strategy to generate point cloud of uniform density.

As described earlier, at this stage we represent the virtual model as a set of triangles of different sizes in a 3D scene. Sampling of points on the triangle surface is performed by random selection of triangle in each iteration and generating a point at random location inside the triangle. If we want to generate a point cloud of size m, we repeat the process of random triangle selection m times, and each time we generate a random point inside the selected triangle. However since the triangle sizes are not uniform, the part of the model having higher number of small triangles will likely to have higher point density than the part having smaller number of large triangles. In order to handle this problem, we use the triangle area as a selection probability while performing the sampling, thus resulting in a list of triangles sampled in a uniform manner.

Next, in order to perform point sampling on each triangle surface, we exploit barycentric coordinates, which are widely used in ray tracing applications of computer graphics. The barycentric coordinates can be used to define the position of any point located inside a triangle. In the general case for a simplex in an affine space, let the vertices are denoted as (v 1 , • • • , v n ). If the following condition holds true for some point P:

(α 1 , • • • , α n )P = α 1 v 1 + • • • + α n v n , (1)
and at least one of the α i 's are not zero, then the coefficients (α

1 , • • • , α n ) are the barycentric coordinates of P with respect to (v 1 , • • • , v n ).
For the case of triangles, any point inside a triangle can be described by a convex combination of its vertices. Let the vertices of a triangle are denoted as v 1 , v 2 , v 3 . Then any random point Q inside the triangle can be expressed as,

Q = 3 i=1 α i v i , (2)
where α 1 , α 2 , α 3 ≥ 0 are three scalars (the barycentric coordinates of the point Q) such that the following condition holds true:

α 1 + α 2 + α 3 = 1. (3)
For all points inside the triangle, there exists a unique triple (α 1 , α 2 , α 3 ) such that the condition of Equation 2 holds true. Hence, due to the property in Equation 3, any point can be defined by any two of the barycentric coordinates.

Using barycentric coordinates, random points can be generated easily on the triangle surfaces. For each randomly selected triangle, we generate two random numbers α 1 and α 2 . Then α 3 is computed as, α 3 = 1-(α 1 +α 2 ) (using Equation 3). This gives us a random point Q inside the triangle. By performing this process for m number of times, a uniformly spaced point cloud of size m is obtained on the surface of the virtual model.

Insideness Testing

Although the point sampling is performed on the surface of the primitives, there can be points remaining inside the volumetric primitives (cylinder, frustum, and sphere). As shown in Figure 1, for the case of a simple branching structure approximated by cylinders and frustums, part of a branch remains inside the other. This is physically inconsistent, and we need to filter out the points inside any volumetric primitive. Hence for each point, we perform a geometric testing for different types of primitives as follows.

Cylinder Primitive Let an arbitraily oriented cylinder has radius r and two endpoints ( -→ p i , -→ p j ) on the central axis. Given a random point -→ q , we aim at determining if the point is inside the cylinder or not (left of Figure 2). First, we test if -→ q lies between the planes of the two circular facets of the cylinder. This is determined by the following two conditions:

( - → q -- → p i ) • ( - → p j -- → p i ) ≥ 0, ( - → q -- → p j ) • ( - → p j -- → p i ) ≤ 0. ( 4
)
If the above conditions are true, we compute the perpendicular distance of the point -→ q from the central axis of the cylinder. If the distance is less than the radius of the cylinder, then the point lies inside. Mathematically the condition is defined as,

|( - → q -- → p i ) • ( - → p j -- → p i )| |( - → p j -- → p i )| ≤ r. ( 5
)
Fig. 2. Insideness testing for the case of cylinder and frustum.

Frustum Primitive Let an arbitraily oriented frustum has base radius r i and top radius r j . The two endpoints on the central axis are defined as ( -→ p i , -→ p j ). Given a random point -→ q , we aim at determining if the point is inside the frustum or not (right of Figure 2).

Similar to the case of cylinder, first we check if the point lies between the planes of the two circular facets of the frustum using Equation 4. Then we project the point -→ q on the line joining -→ p i and -→ p j (the central axis of the frustum) as,

- → q = - → p i + ( - → q -- → p i ) • ( - → p j -- → p i ) ( - → p j -- → p i ) • ( - → p j -- → p i ) * ( - → p j -- → p i ). ( 6
)
The distance between -→ q and -→ q is computed as, d = dist( -→ q , -→ q ). Next we compute the frustum radius r q at the point -→ q . This can be computed as,

r q = dist( - → p i , - → q ) dist( - → p i , - → p j ) * r i + dist( - → p j , - → q ) dist( - → p i , - → p j ) * r j . ( 7
)
If the distance d < r q , then the point lies inside the frustum.

Sphere Primitive Let a sphere centered at point -→ p = (x c , y c , z c ) having radius r. A point -→ q = (x i , y i , z i ) is inside the sphere if the following condition holds true:

(x c -x i ) 2 + (y c -y i ) 2 + (z c -z i ) 2 ≤ r. (8)
General Primitives For more general geometric models that represent more complex plant organs, we use general computer graphics techniques like ray casting [START_REF] Horvat | Ray-casting point-in-polyhedron test[END_REF] that make it possible to test whether a particular point lies inside or outside a closed geometric envelop.

The step by step procedure to generate the point cloud using the basic point sampling strategy is outlined in Algorithm 1.

Sensor Specific Distance Effect

The point cloud sampling described so far, is based on the assumption that the point density is uniform in the whole plant model. However in real scanning devices, this is usually not the case. While there are different types of sensor devices having different types of specific traits associated with the generated surface point cloud, one of the most common is the effect of distance from the scanner on the point density. A typical effect of this type yields higher point density near the scanner, and the density decreases as the distance of different parts of the object from the scanner increases. In order to simulate the effect of these type of real scanning devices, we implemented a sampling strategy which can be adapted to different types of cases according to the need of the application.

Instead of random sampling of triangles in every iteration as in the case of naive sampling technique described before, we compute the number of points to be sampled on each triangle based on a distance factor. First we compute the total surface area A T of all the triangles as, A T = i A Ti . Then we compute the point density ρ u per unit area of the mesh as, ρ u = n/A T , where n is the maximum number of points we wish to sample. Basically ρ u gives the point density for the case of uniform sampling: given a random triangle T i , the number of points on the triangle surface can be approximated as (ρ u * A Ti ).

Next we compute the distance d i of each triangle centroid from the camera position, and normalize the value in the range [0, 1]. Now the idea is that, as the distance of a triangle increases from the camera position, the point density ρ u is decreased by some function of d i . For i -th triangle, the number of points to be sampled is computed as,

n i = ρ u A Ti γ(d i ), (9)
where γ(•) is a real valued function that controls the effect of distance on the point sampling. To demonstrate the effect of γ(•), we show an example of a simple branching structure in Figure 3, where we used γ(d i ) = 1 (which is the uniform sampling), and γ(d i ) = d 2 i . The Figure shows how the point density decreases as the distance of the object part increases from the camera position, compared to the case of uniform sampling.

Effect of Noise

Presence of noise is a very common effect of scanning devices. While the type of noise can vary among different types of devices, we simulate a few common effects of noise in our framework. Some examples are shown in Figure 4 on a simple branching structure as in Figure 3. On the left of the figure, we show the effect of uniform noise, where each point is generated by a Gaussian having the surface position of the point as mean and a standard deviation of σ = 0.06. We also check the insideness testing while applying the noise, so that the noisy points do not remain inside any branch. In the middle of the figure, the noise is modelled with the effect of distance from the camera, as described in the previous section.

Points closer to the scanner has almost zero noise, while the noise increases as the distance of the point increases from the position of the scanner. Finally at the right of the figure, the distance effect is shown considering the decrease of point density as described before.

Results

We have implemented the model in L-Py framework [START_REF] Boudon | L-py: an lsystem simulation framework for modeling plant architecture development based on a dynamic language[END_REF]. However, other types of frameworks (e.g. L+C [START_REF] Karwowski | Design and implementation of the l+c modeling language[END_REF]) can also be used. In the left of Figure 5, we show the synthetic model of an Arabidopsis plant that is generated by parametric Lsystem using stochastic parameters to compute different instances of the plant architecture. Using the stochasticity in the model, it is possible to generate a large variety of plants along with their labelled point cloud data. In the same figure we show some of the varieties of the generated point cloud models of the Arabidopsis plant. In all examples, we have used 3 labels (flowers, internodes and leaves) and about 100k points, although the number of labels and points can be controlled according to the need of the application. In Figure 6, we show more detailed result of the Arabidopsis plant. We also show some results on different types of artificial plants in Figure 7. This demonstrates the portability of the framework which allows us to use it in any procedural plant model. In the last example (rightmost in the same figure), we demonstrate that the number of labels can be controlled according to the need of the application. In this case, we have used different labels for each flower of the Lilac plant (denoted by different colors).

Discussion

There are many ways the framework can be improved and adapted to different types of applications. For example, incorporating back face culling in the framework can be useful for applications where occlusion needs to be modelled. While we have implemented a basic sampling strategy to mimic the sensor specific distance effect, there can be several other types of effects in real world scanning devices, that yield a research avenue by themselves. Incorporating these effects in the point sampler will make the system more versatile and robust, and increase its application range. As we have simulated the effect of distance on the point density by considering the distance of the triangle centroids from the scanner, the point density remains same in all areas of one triangle. For large triangles where one part of the triangle is closer to the scanner than the other part, ideally the point density should vary accordingly inside the same triangle. However in the current framework, we have not modelled this effect.

In order to model the sensor specific distance effect, the desired number of points should be significantly higher than the total number of triangles in the model. In the current implementation, we estimate the number of points per triangle based on its area, without taking into consideration the number of sampled points on the neighboring triangles. While the number of estimated points per triangle can be a fractional number according to Equation 9, we need to round off the value to a nearest integer. Hence the part(s) of the model having large number of tiny triangles will contain at least one point inside each triangle, even when the number of estimated point for the triangle is less than 1, resulting in higher point density in these areas. Ideally the sampling should be performed based on the point density of the neighboring triangles, although this is not considered in the framework.

The current implementation is not optimized, and thus is slow in terms of computation time. For example, in order to generate 100k points in the Arabidopsis model, it takes about 1.5 hours to run the basic sampler. While we have decomposed each shape into triangles to perform point sampling, the points can also be sampled directly on the surface of the geometrical shape. However, we use the triangle representation since it is simple and straightforward.

Conclusion

In this paper we have proposed a generalized approach to generate labeled 3D point cloud data from procedural plant models. As a preliminary work, we have presented results on a simulated Arabidopsis plant along with some artificial plant models. The point sampler is portable enough to use in different types of plant models, and corresponding point cloud data can be generated. This can be beneficial to various types of plant phenotyping research problems in performing quantitative evaluation of point cloud based algorithms, as well as in training the machine learning models. It will be worth investigating how the synthetic models perform in training deep neural networks. In future, we plan to model different types of plant species and more complex plants in the framework. While currently there is no such method available to generate synthetic labeled point cloud data, our approach is the first step towards this in promoting 3D point cloud based plant phenotyping research. Randomly select a triangle T i j ∈ T with area probability {v j 1 , v j 2 , v j 3 } ← vertex(T i j ) // Extract triangle vertices α1 ← rand(0, 1), α2 ← rand(0, 1), α3 ← rand(0, 1) α1 ← α1/ 3 i=1 αi, α2 ← α2/

Fig. 1 .

 1 Fig. 1. In case of volumetric primitives, joining of two (or more) branches results in one primitive to retain inside another primitive.

Fig. 3 .

 3 Fig. 3. Simulation of the effect of distance from the scanner position. From left to right: the original model along with camera position, using γ(di) = 1 (uniform point sampling), and γ(di) = d 2 i .

Fig. 4 .

 4 Fig. 4. Simulation of noise effect. Left: uniform distribution of noise, Middle: Effect of noise on the distance from the scanner, Right: Effect of noise on the distance along with point density factor.

Fig. 5 .

 5 Fig. 5. From left to right: a virtual model of Arabidopsis plant, followed by labelled point cloud representation of the plant model in different growth stages.

Fig. 6 .

 6 Fig. 6. An example of a point cloud model generated by the approach using 3 labels and about 100k points. Different parts are shown in details by the zoom-in.

Fig. 7 .

 7 Fig. 7. Some examples of artificial plant models (top row), and the corresponding labelled point cloud data (bottom row). In the first two examples, 2 and 3 labels are used, respectively. In the last example, each flower is assigned a different label denoted by different colors.

Algorithm 1 :

 1 Basic Point Sampling Algorithm Input: Procedural Plant Model Output: Surface point cloud of the model P = {p1, • • • , pn}, label list begin P ← {Φ}, ← {Φ} Assign label i, i ∈ {1, • • • , r} to r organs in the procedural model Extract all primitives with organ labels from the model {C i 1 , • • • , C i m } Extract the set of geometry parameters of the primitives {λ1, • • • , λm} Tessellate each primitive into set of triangles T = {{T i 1 }, • • • , {T i m }} Create global list of triangle vertices along with their corresponding labels while Total number of Points ! = n do

  Cylinder then insidenessFlag ← isInsideCylinder(λ k , p)// using Eqn.4,5 else if C i k is a Frustum then insidenessFlag ← isInsideFrsutum(λ k , p)// using Eqn.4,6,7 else if C i k is a Sphere then insidenessFlag ← isInsideSphere(λ k , p)// using Eqn.8

	p ← α1v j 1 + α2v j 2 + α3v j 3	3 i=1 αi, α3 ← α3/ 3 i=1 αi
	insidenessFlag ← F alse
	// Loop over all primitives
	for k ← 1 to m do
	if C i k is a if insidenessFlag = False then
	P ← P ∪ p
	← ∪ i

Acknowledgment

This work is supported by Robotics for Microfarms (ROMI) European project.


Damienne Provitolo
email: provitol@geoazur.unice.fr

René Lozi
email: rene.lozi@univ-cotedazur.fr

E Tric

Topological analysis of a weighted human behaviour model coupled on a streets and places network in the context of urban terrorist attack

Keywords: human behaviour, panic, terrorism, street and place network, urban context, mathematical modelling, Coupled Weighted PCR model

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the context of disasters, and in order to better protect the population, one of the major challenges today is to better understand and anticipate both individual and collective human behaviour, and the dynamics of the displacements associated with these behaviours. Indeed, the impact of a dangerous phenomenon is particularly determined by the behaviour of the affected population. These reactions allow anyone to ensure their own safety and that of their family. This is especially true in the case of sudden and unpredictable events, such as terrorist attacks. These events require immediate reactions for self-protection and self-evacuation, before the arrival of emergency response services.

In spite of the long history of terrorism, there is currently no uniform definition of this word, because "depending on the political configuration, this one will be terrorist for some; instead he will be hailed freedom fighters for others" [START_REF] Marchesin | Introduction aux relations internationales[END_REF]. However, it is possible to identify different forms of terrorist acts, such as attacks, kidnappings, sabotage, bio-terrorism, or cyber-terrorism. Moreover, terrorist acts can be committed at different levels by individuals, groups or states.

In this article, we focus specifically on the new forms of terrorism, such as the attacks perpetrated by sects or groups which are becoming more and more active. These attacks target places frequented by the public (bus and subway stations, airports, hotels, cafes, etc.), in order to spread panic among the population. The sarin gas attack on the Tokyo subway committed by the Aum sect (1995), the bomb attacks in bus and subway stations (Nigeria, 2014; Belgium, 2016), at airports (Brussels airport, 2016; Atatürk airport, 2016) or near the finish line of the Boston Marathon (2013), and acts perpetrated by Islamic organizations (Kashmir, 2019), are a few examples of a very long list of terrorist attacks. All these events triggered mass movements of flight and collective panic.

Moreover, in order to cause as many casualties as possible and to complicate the response of the security forces, terrorists have developed a new modus operandi by conducting simultaneous attacks in urban areas. In this way, for some coordinated attacks, human losses have been massive, yet constrained in time and space. This was the case for September 11, 2001, in New York and near Washington D.C., where a series of four coordinated terrorist attacks claimed by the al-Qaeda terrorist group killed nearly 3,000 victims and injured over six thousand. India also endured terrifying terrorist attacks in November 2008 when an Islamic terrorist organization based abroad carried out a series of 12 coordinated shooting and bombing attacks lasting four days across Mumbai. This modus operandi also makes it possible to prolong the attacks to maximise media coverage and to increase the feeling of terror and panic, as was seen during the Paris and Saint-Denis (France) terrorist attacks of November 2015.

In response to tightened police controls to reduce bombing attacks, new types of ad hoc weapons are now employed, like lorries or cars simply driven directly through a crowd massed for a cultural event or simply strolling peacefully. These attacks are most often led by "lone wolves", who are difficult to detect beforehand and equally difficult to locate because they are "nested" among civilians. On July 14, 2016, 86 people who attended the Bastille Day fireworks on the "Promenade des Anglais" in Nice, France, were killed; 458 more were wounded. Since this day the frequency of such kinds of terrorist attacks has grown rapidly in western countries: in Berlin, Germany, December 19, 2016; in London, England, March 22, 2017; in Stockholm, Sweden, April 7, 2017; and in Barcelona, Spain, August 16, 2017. It is not only cars that have been used in such attacks, but true weapons like military guns (Las Vegas, USA, October 2, 2017) or simple tools like knives (Marseille, France, October 1, 2017), hammers, etc.

Nowadays, in order to protect inhabitants, authorities have developed diversified risk reduction strategies. These are complementary, with some acting directly on the threat, and others targeting the vulnerability of potential victims. In the field of counter-terrorism, the most widespread strategies aim to fight against the threat. Thus, in order to avoid new attacks, actions aiming to dismantle terrorist networks are carried out by the police and intelligence services. Attention is therefore focused on dangerous groups, but it is difficult to detect and neutralise all terrorist threats, especially those perpetrated by such "lone wolves". Thus, prevention policies also focus on reducing the vulnerability of the population. It is then a question of hindering the terrorist action itself by deploying security measures as close as possible to the potential targets (festivals, concerts, college, train and metro stations, etc.). In many countries, like in India and France, more law enforcement officers and even the army have been deployed for such measures; however this can be done only for a limited period of time.

More recently, in France, other actions, such as the distribution of leaflets and information documents have been deployed, to increase awareness among the population about the responses that could save lives, for example to: flee where possible, barricade the entrance, hide behind a solid obstacle, turn off the phone, etc. This effort of dissemination of the best practices in the face of an attack is already a step forward but it remains incomplete.

We must also look at the reality of human behaviour (i.e. what people actually do during a terrorist attack). Their reactions simultaneously depend on their own emotions, the culture of risk and the environmental context (a closed environment like an auditorium or theatre, compared to an open geographical environment, such as a public square, place or networks of streets). When such an attack happens, the topography of the area is very important. The dynamics of human reactions and the associated displacements are guided by the space and the alternatives that it offers, especially in terms of evacuation, flight, and accessibility to refuge areas.

Consider, for example, the terrible Jallianwala Bagh massacre, which took place on April 13, 1919 in Amritsar, Punjab, when a crowd of non-violent protesters were fired upon by troops of the British Indian Army.

The Jallianwalla Bagh is a public garden with an area of 28,000 m 2 , walled on all sides with five entrances (Figs. 5.1 and 5.2). The largest entry point was blocked by a tank and the main exit was locked. The troops fired on the crowd, directing their bullets largely towards the few open gates through which people were trying to flee (Fig. 5.3). In an urban context, it is also important to study the impact of the width of alleys and the shelter capacity of public squares or places so that the population can ensure its own selfevacuation (Fig. 5.4). Of course, it is difficult to artificially reproduce a disaster, which would otherwise allow us to observe the diversity of human reactions that could occur, to follow the spatio-temporal dynamics and to analyse the impact of territorial configurations on these dynamics. To overcome these limits, it is possible to develop mathematical models from which evolution scenarios are simulated by varying parameters or initial conditions. The computer thus becomes the virtual laboratory and the simulation is understood to be an experiment on a model, a digital experience [START_REF] Colonna | Expériences virtuelles et virtualités expérimentales[END_REF].

The Com2SiCa research team 1 proposed the Panic-Control-Reflex (PCR) model [START_REF] Provitolo | Les comportements humains en situation de catastrophe: de l'observation à la modélisation conceptuelle et mathématique[END_REF], which is a model that simulates the possible human behaviour that can occur during sudden onset and unpredictable disasters, such as a tsunami, earthquake, or technological disaster. This model is formalized by a system of ordinary differential equations to describe behavioural dynamics over time [START_REF] Verdière | Mathematical Modeling of Human Behaviours During Catastrophic Events[END_REF][START_REF] Verdière | Understanding and Simulation of Human Behaviours in Areas Affected by Disasters: From the Observation to the Conception of a Mathematical Model[END_REF][START_REF] Cantin | Mathematical Modeling of Human Behaviours During Catastrophic Events: Stability and Bifurcations[END_REF].

In this article we propose an extension of the PCR model in order to take into account the influence of spatial configuration in the mathematical modelling of the dynamics of human reactions in the face of traumatic situations, such as terrorist attacks. This dynamic and the associated displacements are indeed guided by the territorial configurations (networks of streets and places) and the alternatives that they offer in terms of evacuation, flight, and accessibility to shelters. We call this extended model the Coupled Weighted PCR (CWPCR).

In Section 5.2, we will present the Panic-Control-Reflex model in its graphical and mathematical formalism, as published by the authors cited above. Then, in Section 5.3, this model will be improved in order to take into account the role of spatial configurations on behavioural dynamics. The CWPCR model incorporates the pressure and counter pressure of the crowd in each place via a bottleneck effect, which is induced by the narrowness and the length of the streets and the size of places. In Section 5.4 we will consider an oriented network with three nodes representing three places or public squares of different sizes, linked by narrow streets. This part is therefore devoted to the analysis of the impact of the parameters (the size of the places, the width of the streets) on the evacuation of the population in the face of a terrorist attack, by means of numerical simulations. The numerical results highlight that, depending on their respective size, intermediate places modulate the dynamics and the speed of flow of the crowd. In this sense, they become strategic places both for the planners who must think about the organization of the area to host public events and festivals, and also for the terrorists who can use these strategic places to multiply the effect of their harmful actions by trapping the flight movements between two areas of action. This model is used with a user-friendly graphical representation, which allows planners to accurately consider where to organize host public events in a specific territorial context. Finally, in Section 5.5, a brief conclusion will be drawn.

The PCR system: an a-spatial model for analysing the dynamics of human behaviour during a disaster

Neuroscientific background of the PCR system

The PCR model [START_REF] Provitolo | Les comportements humains en situation de catastrophe: de l'observation à la modélisation conceptuelle et mathématique[END_REF][START_REF] Verdière | Mathematical Modeling of Human Behaviours During Catastrophic Events[END_REF][START_REF] Verdière | Understanding and Simulation of Human Behaviours in Areas Affected by Disasters: From the Observation to the Conception of a Mathematical Model[END_REF][START_REF] Cantin | Mathematical Modeling of Human Behaviours During Catastrophic Events: Stability and Bifurcations[END_REF] is a simulation model of the dynamics of collective human behaviour during a disaster (Fig. 5.5). It has been developed from the SIR based-models, which are compartmental models that are widely used in epidemiology [START_REF] Murray | Mathematical Biology I: An Introduction[END_REF]. In these models, the population can be decomposed into several subpopulations, each of which corresponds to a compartment. The PCR model focuses on: i. different human behavioural states, namely daily behaviour before a disaster occurs, as well as reflex, panic and controlled behaviours that are observed during a disaster, ii. transitional processes from one reaction to another. Indeed, neuroscience research shows that in disaster situations, humans are rarely stuck in one type of behaviour. The population switches between different behavioural states, some of which are the result of instinctive reactions [START_REF] Laborit | La légende des comportements[END_REF], others of reasoned reactions [START_REF] George | Premières impressions[END_REF].

To take these behavioural sequences into account, the PCR model formalizes human reactions as a chain of behaviours that appear in a certain order. It distinguishes q(t), the daily behaviour before, and b(t), after the disaster; r(t) and p(t), the uncontrolled emotional behaviours which are managed by the reptilian zone of the brain; and c(t), the reasoned behaviours that are controlled by the prefrontal cortex [START_REF] Noto | Médecine de catastrophe[END_REF][START_REF] Soussignan | Un monde d'émotions[END_REF][START_REF] Brosch | Les effets cognitifs des émotions[END_REF]. This is represented in Fig. 5.5.

As the brain switches from one behavioural state to another, in the context of terrorist acts and therefore in a situation of sudden and unforeseen threat, the whole impacted population first adopts a behavioural reaction, Reflex r(t), under the influence of surprise and the suddenness of the event, before transiting to the Panic reflex behaviour p(t), or Controlled behaviour c(t).

Reflex r(t) and Panic p(t) behaviour are instinctive, automatic reactions, allowing one to react extremely quickly to the threat, either by being stunned and paralysed r(t), or by fleeing as quickly as possible due to the panic fear p(t). In a context of a dense crowd, the context sought by terrorist groups, panic escape behaviour can worsen the vulnerability of the population because of the risk of crushing and suffocation [START_REF] Crocq | Les paniques collectives[END_REF].

Controlled behaviour c(t) concerns reasoned and self-control reactions. They can take different forms during a catastrophe, for example, in the form of evacuation, leak, containment, sheltering, search for help, mutual aid or, on the contrary, looting, etc. Despite their diversity, the PCR model aggregates all of these controlled behaviours.

During the event, the switches from one behavioural state to another are caused by transitional dynamics due to: i. Causal relationships (B1, B2, C1, C2). Once the population is in the reflex behaviour state, a part of it can evolve towards controlled behaviours at the rate B1, while another part transitions towards panic behaviours at the rate B2.

Likewise, a part of the panicked population may switch to controlled behaviour at the rate C1. According to the evolution of the situation, individuals who have adopted a controlled behaviour may switch back to panic behaviour at the rate C2. This process can be iterated many times. ii. Processes of imitation and contagion, which are well-known in crowd psychology and have been termed "emotional contagion" [START_REF] Hatfield | Emotional Contagion[END_REF]. The imitation processes are modelled identically to epidemiological propagation [START_REF] Provitolo | Un exemple d'effets de dominos : la panique dans les catastrophes urbaines[END_REF]. The imitation is valid in both directions and is modelled by the function F(r,c) Eq.( 5.3) for emotional contagion between reflex and controlled behaviour (using the damping coefficients α1 and α2), by the function G(r,p) Eq.( 5.4) for emotional contagion between reflex and panic behaviour (δ1 and δ2), and by the function H(c,p) Eq.( 5.5) for emotional contagion between controlled and panicked behaviour (μ1 and μ2). iii. Domino effects, which illustrate a succession of events (s1 and s2) corresponding, for example, to a new attack in urban area or to a closed door during an evacuation. In the PCR model, the parameters s1 and s2 are either constant or built in a periodic form.

The triggering of the threat is represented by a forcing function (t), the form of which may vary according to the specificities of the danger (event with fast or slow kinetics, expected or not).

Equations of the PCR system

In [START_REF] Verdière | Mathematical Modeling of Human Behaviours During Catastrophic Events[END_REF][START_REF] Verdière | Understanding and Simulation of Human Behaviours in Areas Affected by Disasters: From the Observation to the Conception of a Mathematical Model[END_REF][START_REF] Cantin | Mathematical Modeling of Human Behaviours During Catastrophic Events: Stability and Bifurcations[END_REF] the authors introduce the Panic-Control-Reflex model (PCR) by using a system of five Ordinary Differential Equations (ODE), which describe the human behaviours in one specific place during a catastrophic event. We include these ODE Eqs. (5.1)-(5.6) in order to describe the modifications we introduce to them in the next Section, where we seek to model such behaviours when a network of places and streets or stairs linking those places in a town is considered. In Section 5.3 we will identify this network to a mathematical graph, the places are called vertices or nodes, and the streets and stairs are the edges linking those nodes. and  given by ( )

( , ) ( , ) T i t X t X  =  , 1,...,5 i = ,
where the functions i  are defined by

1 1 2 1 2 2 1 1 2 1 () ( , ) ( ) ( ) 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( , ) ( ) ( )(1 ( )) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( m rt t X t q t B B

r t s t c t s t p t r F r t c t r t c t G r t p t r t p t t X t c t b t B r t C p t C c t s t c t F r t c t r t c t H c t p t c t p

    = - - + + +   ++  = - - + + - - -+ 3 2 1 2 2 4 5 ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) () ( , ) ( ) ( ) 1 ( , ) ( ) ( )(1 ( )) m t t X B r t C p t C c t s t p t G r t p t r t p t H c t p t c t p t rt t X t q t r t X t c t b t            = - + - -  -     = - -      = -  (5.2)
and the variables ( ), ( ), ( ), ( ), ( ) r t c t p t q t b t denote, respectively, the densities of people being in a reflex, control, panic, daily or back to daily behaviour [START_REF] Cantin | Mathematical Modeling of Human Behaviours During Catastrophic Events: Stability and Bifurcations[END_REF].

The parameters involved in (5.2) are real positive coefficients previously defined in Sec. 5.2.1:

0 m r  (reflex behaviour maximum value); 0, 0, 1, 2 ii B C i   = ; 0, 0, 0, 1, 2 i i i i       = ; 0, 1, 2 i si =.
The imitation functions , FG and H are real valued functions defined on ℝ × ℝ by

1 1 2 2 ( ) ( ) ( ( ), ( )) ( ) ( ) r t c t F r t c t f f c t r t       = - +     ++     (5.3) 1 1 2 2 ( ) ( ) ( ( ), ( )) ( ) ( ) r t p t G r t p t g g p t r t       = - +     ++     (5.4) 1 1 2 2 ( ) ( ) ( ( ), ( )) ( ) ( ) c t p t H c t p t h h p t c t       =-    ++     , ( 5.5)
where  is a positive number and ,,

i i i f g h for 1, 2 i =
are real valued functions defined on ℝ. They have a decreasing shape indicating that the behaviour imitation is symmetric. Moreover they are normalized,

0 ≤ 𝑓 𝑖 (𝑢) ≤ 1, 0 ≤ 𝑔 𝑖 (𝑢) ≤ 1, 0 ≤ ℎ 𝑖 (𝑢) ≤ 1, ∀𝑢 ∈ ℝ, 𝑖 = 1,2 (5.6)
Because this model does not take the mortality rate into account, the population remains constant and, in one node, can be normalized to 1. Therefore Eq. (5.1) is considered when time is proceeding from an initial time 0 0 t  , with initial condition 0 0 0 0 0 0 0 0 0 0 ( ( ), ( ), ( ), ( ), ( )) ( , , , , ) r t c t p t q t b t r c p q b = (5.7) that satisfies the following properties 0 0 0 0 0 ( ) 0, ( ) 0, ( ) 0, ( ) 0, ( ) 0

r t c t p t q t b t     
(5.8)

0 0 0 0 0 ( ) ( ) ( ) ( ) ( ) 1 r t c t p t q t b t + + + + =
(5.9) Equation (5.9) remains true throughout the process because the sum of the five equations (5.2) is null, therefore

0 ( ) 1 ( ( ) ( ) ( ) ( )), b t r t c t p t q t t t = - + + +   (5.10)
which implies that Eq. ( 5.2) can be reduced to the system of only four ODE

1 2 1 2 1 1 2 1 () ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ) ( ) ( )(1 ( )) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) m rt r t t q t B B

r t s t c t s t p t r F r t c t r t c t G r t p t r t p t c t t c t b t B r t C p t C c t s t c t F r t c t r t c t H c t p t c t p t p

   = - - + + +   ++ = - - + + - - -+ 2 1 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) () ( ) ( ) ( ) 1 m t B

r t C p t C c t s t p t G r t p t r t p t H c t p t c t p t rt

q t t q t r          = - + - -  -    = - -     (5.11)
The initial condition of Eq.(5.11) corresponding to 0 0 0 0 0 ( , , , , ) r c p q b for Eq. (5.2) becomes simply 0 0 0 0 ( , , , ) r c p q

(5.12)

Transitional dynamics

Both forcing functions,  and  , respectively model the beginning of the disaster and the return to a quiescent daily behaviour. Their shape can be adapted to various scenarios. When t is sufficiently large they satisfy ( ) ( ) 1 tt  == . In catastrophic situations it is considered that  is a stiff function, ranging from 0 to 1 in a very brief interval of time [START_REF] Verdière | Mathematical Modeling of Human Behaviours During Catastrophic Events[END_REF][START_REF] Verdière | Understanding and Simulation of Human Behaviours in Areas Affected by Disasters: From the Observation to the Conception of a Mathematical Model[END_REF][START_REF] Cantin | Mathematical Modeling of Human Behaviours During Catastrophic Events: Stability and Bifurcations[END_REF] because if we consider a bomb attack, all the crowd that is near the explosion passes from daily to reflex behaviour in an instant, and it takes a very long time for people to return to their normal state. Hence, during the transitional dynamics, the population with daily behaviour collapses and there is not yet a population that is back to daily behaviour (i.e. ( ) ( ) 0 q t b t == ).

    -   =          , (5.14) 2 2.5 cos 2 , 0 2.5 () 10 1, 2.5 x if x t if x    -     =      , ( 5
System (5.11) is reduced to

1 2 1 2 1 1 2 1 2 1 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ((

r t B B r t s t c t s t p t F r t c t r t c t G r t p t r t p t c t B r t C p t C c t s t c t F r t c t r t c t H c t p t c t p t p t B r t C p t C c t s t p t G r t

= -+ + + + + = + - - - + = - + - - ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) . p t r t p t H c t p t c t p t        -
(5.17)

Mathematical Weighted and Coupled PCR system

The Weighted PCR system

Equations (1-4) model the human behaviours, in one specific place, during a catastrophic event. As explained in the first part of this article, it is based on neuroscience studies.

However just before a dramatic event like a terror attack in a city, the crowd is generally spread across several places, public squares and streets. In the aftermath of the initial shock, people are rushing through the streets to reach what they think will be more secure places (Fig. 5.7). linking these nodes. They are oriented, because the flow from one place toward another is not symmetric. Our aim is to model the motion of the crowd through such edges. To achieve this, we must to introduce some "geographical" particularities of the city, like the size of places, the narrowness of streets, and the number of people initially present in each place. This is why we need to we call ( ), ( ), ( ), ( ), ( )

k k k k k
r t c t p t q t b t the number of people being in reflex, control, panic, daily and back to daily behaviour and () k Vtthe total number of people present at this node

( ) ( ) ( ) ( ) ( ) ( ) k k k k k k V t r t c t p t q t b t = + + + + (5.18)
Since we consider that the nodes are not identical (as is seen in an actual city) they do not generally contain, at each moment, the same number of people. Moreover, this number is varying with time when the crowd is moving through the streets. That is why it is more convenient to consider that the five variables of the WPCR model represent actual numbers of people, rather than densities, as in PCR model. Densities can be used only when there is no motion at all between places, and the population in each place, is the same. Of course, this number of people can be transformed locally to density whenever it is necessary.

Secondly, in order to more precisely model the characteristics of the city, we introduce k W , the maximum capacity of the number of people who can be present in each node k N (i.e. due to the size of the corresponding place). Each maximum capacity is a constant. At every time one must have ()

kk V t W  (5.19)
The Weighted Panic-Control-Reflex model (WPCR) is then defined on each node k by

T k k k k k k X r c p q b R =
and  given by ( )

( , ) ( , ) T k i k t X t X  =  , 1, ..., 5 i =
, where the functions i  are defined by ( )

1 1 2 1 2 2 1 1 2 1 ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( , ) ( ) ( ) ( )( ( )) ( ) ( ) ( ) ( ) ( ) ( ( ), k k k k k k k k k k k k k k k k k k k k k k k k k kk t X

r t t q t W r t B B r t s t c t s t p t F r t c t r t c t G r t p t r t p t t X c t t c t W b t B r t C p t C c t s t c t F r t c

   = = - - + + + + +  = = - - + + - - - ( ) 3 2 1 2 2 4 5 ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ( ), ( )) ( ) ( ) ( ( ), ( )) ( ) ( ) ( , ) ( ) ( ) ( ) ( ) ( , ) ( ) ( ) ( )( k k k k k k k k k k k k k k k k k k k k k k k k k k k k t

r t c t H c t p t c t p t t X p t B r t C p t C c t s t p t G r t p t r t p t H c t p t c t p t t X q t t q t W r t t X b t t c t W  

+  = = - + - - -  = = - -  = = ( )) kk bt           -  (5.21) the initial condition satisfies 0 0 0 0 0 0 ,0 ,0 ,0 ,0 ,0 ,0 ( ) ( ) ( ) ( ) ( ) ( ) k k k k k k k k k k k k k r t c t p t q t b t V t r c p q b V W + + + + = = + + + + =  (5.22)
We suppose that the characteristic parameters 1 2 1 2 , , , , B B C C of each node have the same value, because they depend on cultural and psychological factors specific to each individual rather than to spatial configurations and crowd context. Thus, all the parameters and the functions are the same as those defined in Sect. 5.2.2., which is why we use function i  instead of function , ki  .

Flows and bottleneck coupling

As previously defined, nodes are linked by edges. We now aim to model the motion of the crowd through such edges (i.e. streets, stairs, escalators, doors, etc.). For the sake of simplicity, we present the coupling on a simplified oriented network with only two nodes . We suppose that during the short interval of time when people are travelling inside one edge, they remain in the same behavioural class (Fig. 5.8.). In this figure, the "geographical" edge 12 () NN → (i.e. the street linking node 1 to node 2) is split in three "behavioural edges" meaning that on the same street people in reflex, panic or controlled behaviour are escaping from node 1 to node 2, therefore people in each particular behaviour in node 1 are meeting people in the same behaviour in node 2. To continue focusing on the special coupling that we are introducing here, we suppose that imitation mechanisms are not activated (i.e. 0 F H G    , which is equivalent to

0 i i i    = = = for 1, 2 i =
), and, furthermore, there is no domino effect ( 0, 1, 2 i si ==) and we also suppose that we are in the interval   such transitional dynamics are the solution of the system which is the reduction of system (5.21) in the transition interval (as Eq. (5.17) is the reduction of Eq. (5.2)).

12 1 1 2 2 1 2 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) kk k k k k k k k r t B B

r t c t B r t C p t C c t p t B r t C p t C c t

= - +   = + -   = - +  (5.
Note: in the WPCR model, the terms ( ), ( ), ( ), ( ), ( )

k k k k k
r t c t p t q t b t can be considered as flows, because a flow is a quantity of something divided by a unit of time. There are two kinds of flow. In Eqs. (5.21)-(5.23) flows are "behavioural", as they represent the quantities of people changing their behaviour per unit of time. Now we consider also "motion" flows, which are the quantities of people in each behaviour class, moving from one node to another node, per unit of time. Of course, both kinds of flow are combined to produce a global equivalent in the following equations.

In [START_REF] Vanumu | Fundamental diagrams of pedestrian flow characteristics: A review[END_REF] different types of flow situations are considered in pedestrian facilities, such as unidirectional, bidirectional and crossing: it is said that "Unidirectional and bidirectional flow conditions can be commonly observed in corridors, stairs and bottlenecks of pedestrian facilities such as transport terminals, shopping malls, etc. An understanding of the fundamental relationship between flow -speed -density is important in the planning, design and operation of pedestrian facilities. Capturing the realistic behaviour of pedestrians in various pedestrian facilities with different geometric elements such as corridors, bottlenecks, stairs, escalators etc. is essential in order to estimate the flow parameters accurately. The important parameters such as width of the bottleneck slope of the stairs plays a vital role in deciding the capacity of the respective element. The flow density relationship for different geometric elements is important and further analysis like spatial and temporal development of the basic quantities (velocity, density and flow) on different elements like corridors, stairs and bottlenecks should be considered". Many studies on pedestrian flows have been published [START_REF] Daamen | First-Order Pedestrian Traffic Flow[END_REF][START_REF] Dias | Investigating collective escape behaviours in complex situations[END_REF][START_REF] Kretz | Experimental study of pedestrian flow through a bottleneck[END_REF][START_REF] Liao | Experimental study on pedestrian flow through wide bottleneck[END_REF][START_REF] Seyfried | The Fundamental Diagram of Pedestrian Movement Revisited[END_REF][START_REF] Zhang | Empirical study of a unidirectional dense crowd during a real mass event[END_REF][START_REF] Zhou | Pedestrian Evacuation Time Model for Urban Metro Hubs Based on Multiple Video Sequences Data[END_REF]. We consider, in particular, the survey of [START_REF] Daamen | First-Order Pedestrian Traffic Flow[END_REF], in which the graphs of six different experiments showing the relationship between the flow of pedestrians and their densities are displayed. All these graphs show clearly a nonlinear relationship of a logistic type between density and flow (Fig. 5.9). Moreover, Daamen and co-authors developed a first-order traffic flow theory to describe two-dimensional pedestrian flow operations in the case of an oversaturated bottleneck in front of which a large, high-density region has formed (Fig. 5.10). Such a mathematical model also highlights the logistic relationship for any bottleneck width. Fig. 5.10 Relationship between flow and density of pedestrians going through an oversaturated bottleneck, from [START_REF] Daamen | First-Order Pedestrian Traffic Flow[END_REF] We now introduce our hypothesis for the coupled WPCR (CWPCR), based on this type of logistic relationship. Again, for the sake of simplicity, we consider only unidirectional motion of the crowd (i.e. motion on an oriented graph) as in Fig. 5.8, because it is supposed that a terrorist attack occurs in node 1 and that people try to escape from this node towards node 2. As an aside, in forthcoming research we will allow bidirectional motion in more complex network, where bidirectional motion will be simply obtained by adding symmetric terms in Eq. (5.24).

When the crowd is moving from one node to another, its speed and the corresponding motion flow depends on three factors. The first factor reflects the narrowness and the length of the street. More people can go from one place to the next if the street is large, rather than in the case of a narrow street. This topological characteristic will be modelled by a "roughness" coefficient 1,2  . In fact, in the considered coupling, we suppose that people cannot change their behaviour when they move from one node to another (e.g. controlled people remain controlled, panicked people remain panicked, and so on), and thus we use three such roughness coefficients: ,

  

to refer to the population in reflex, panic or controlled situations, respectively. Of course, they can have the same value. The second factor is proportional to the number of people present in the first node. This is equivalent to the pressure of the crowd in the first place. The third factor reflects the counter pressure due to the maximal capacity of the second place conjugated with the number of persons already present there. The combination of pressure and counter pressure gives a bottleneck effect.

Fixed points of the two-node system

The fixed-point research allows us to identify the point of equilibrium towards which the system tends during the transitional period. This equilibrium point highlights the primordial role of the size of node 2 in the context of evacuation dynamics. It is important to note that the mathematically calculated equilibrium point does not necessarily correspond to the equilibrium situation sought by crisis management personnel. Thus, the equilibrium point obtained in situation 2 below, in which a part of the population cannot escape and remains in a dangerous place, does not correspond to a crisis equilibrium situation. The population stranded in the initial place (node 1) remains very vulnerable to the terrorist threat.

The fixed point ( * * * * * *

1 1 1 2 2 2
, , , , , r c p r c p ) of the system toward which the solution Eq. (5.24) converges, is easily computed.

* 1 * 1 * 1 * 2 1 1,0 2,0 * 2 12 2 1,0 2,0 * 2 12 * 1 1 1,0 2,0 2 * 1 12 2 1,0 2,0 2 * 1 12 * 2 * 12 2 12 * 22 2 12 0 0 0 0 ( ) ( ) 0 ( ) ( ) 0 r c p r C V V c CC C V V p CC r C V V W c CC C V V W p CC r CW c CC CW p CC  =  =  =   =    +   =  +    +   =  +     =   +-   =  +    +-  =    +   =    =  +   =  +   1,0 1,0 1,0 2,0 2,0 2,0 1,0 2,0 2 1,0 1,0 1,0 2,0 2,0 2,0 1,0 2,0 2 if , if , r c p r c p V V W r c p r c p V V W + + + + + = +  + + + + + = +       (5.28)
The values of this fixed point mean that: Situation 1: if the number of people initially staying in both nodes is less than the capacity of node 2 (i.e. 1,0 2,0 2

V V W +
), after a while, node 1 becomes empty and all the crowd has sought refuge in node 2.

Alternatively, in Situation 2: if this number is greater than the capacity (i.e. Of course, across the world there are different cultures, which lead to different behaviours. These behaviours can be modelled varying parameters.

For example, if populations are not made aware of major risks and not prepare for them, it causes a panic reaction (B1 < B2), this behaviour is then regulated by the ratio 1 2

C C

. The higher this ratio, the more the population remains or transits in the controlled state.

In both simulated situations, we choose a set of parameters that highlight a weak risk culture, while favouring the return to a controlled behaviour: instead, the values, 12 0.2, 0.4 BB == , mean that there is a weak risk culture and 12 0.3, 0.1 CC == , mean that the Examining both Figs. 5.12 and 5.13, it appears that one of the main trends of the CWPCR model is to take into account the role of the spatial configuration on both behavioural dynamics and on the crowd' ability to escape from a dangerous place towards a place of shelter, when such a place offers sufficient room for the crowd. In effect, in the first situation, everyone can leave place 1, as we can see on Fig. 5.12 where V1(t) vanishes for t > 15 mn, and the entire population escapes to place 2 (V2(t) = 800 for t > 15 min). However, in the second situation, a part of the population remains stranded in place 1, because W2 = 500, which is less than the total number of the crowd equal to 800. Therefore, we can see in Fig.

Scaling the parameters

Fixing the value of all the parameters in PCR, WPCR or CWPCR is a very complicated task, which has not yet been done on an experimental basis. In the framework of the Com2SiCa project, an experimental protocol is under scrutiny, in order to achieve these results in the near future. However, it is important to note that only relative values between parameters are important, because there is a relationship between the unit of time and the unit used for the parameters. In other words, considering all the parameters of the adimensional Eq. (5.29), (i.e. , , , ,

i i c i j p i j r i j B C i j      = = , , 0, 0, , , , 1, 2, 2, 3
) integrated with respect to the variable time t , it is nearly equivalent to consider such parameters multiplied by the same constant  and integrated using the time variable t   = (e.g. can be considered in seconds, minutes, or hours). There is not, strictly speaking, equivalence because such parameters are linearly used in the WPCR model, however in the CWPCR model the coupling is nonlinear and a slight distortion intervenes during a transient short period for some variables.

Speed of convergence towards the fixed points

To stay closer to reality, we further assume that people in reflex situation are stunned and paralysed. They cannot rush from one node to another (Fig. . This value means that, initially, all the people are in node 1 and both nodes 2 and 3 are empty, but after a certain period of time, both nodes 1 and 2 are empty and everyone has reached node 3. One can see the flow of people through the two edges 1 [START_REF] Vanumu | Fundamental diagrams of pedestrian flow characteristics: A review[END_REF] Convergence towards the fixed point ( 0, 0, 0, 0, 0, 0, 0, 15000, 5000 )

Influence of the intermediate place capacity on the evacuation dynamics

Intermediate places play a central role in the fluidity or, to the contrary, the congestion of movement between a dangerous place and a shelter place. This can be shown in the following numerical experiments: with the same parameter values, except for the size of node 2, we analyse the speed at which the people are emptying the place of the terrorist (node 1). We consider the following values of (Fig. 5.17). In the case W2 = 1000 (black curves) the flight of the entire population from place 1 to place 2 (which can be a small square) and then to place 3 is very fast, it lasts less than 10 minutes (Figs 5.17c, f, i), because place 3 can foster the entire population. There is a massive influx of panicked people (Fig. 5.17d), which is greater than the controlled one (Fig. 5.17e), into place 2, which empties very quickly as the majority of controlled people reach the safe shelter (Fig. 5.17h). However, it can be noted that a significant number of the panicked population remains in the refuge place, and this number is only slowly decreased (see bump (Fig. 5.17g)). This is explained by the fact that the flight dynamics are not hindered by obstacles or bottlenecks, and the fleeing populations have not enough time to change their behavioural state.

On the other hand, if W2 = 50 (red curves), the evacuation of the total population from place 1 to place 2 and then place 3 is much slower (about 25 minutes instead of less than 10 minutes, (Figs 5.17c, f, i)). The panicked population has time to calm down in place 1, because there is no new attack (it has been assumed that there is no domino effect, s1 = s2 = It can be highlighted from this first analysis that the faster the speed of change of location (hence decreasing the vulnerability of populations), the faster this speed leads to significant flows of panic in both places 2 and 3. There is a paradox here: the fast self-safety movement of populations leads to situations of collective panic that are more difficult to manage. This fact must be taken into account by emergency services and emergency physicians.

Although the deaths that occur are not included in this CWPCR model version, one can imagine that the escape of panicked populations would induce more victims.

As said before, the simulation results shed light on the importance of the size of the intermediate places and their role in the fluidity or, to the contrary, on the congestion of movements between a dangerous place and a shelter place.

Depending on their respective size, intermediate places will modulate the dynamics and the speed of flow of the crowds. In this sense, they become strategic places both for the planners, who must think about the organization of the area to host public events, and also for the terrorists who can use these strategic places to multiply the effect of their harmful actions by trapping the flight movements between two areas of action. This can be summarized in both Figs. 5.18 and 5.19, where the duration of evacuation time for 80% of people from place 1 is displayed according two variables: the capacity of place 2 versus that of place 3. The arrow goes from long durations (in warm colours) to short durations (in cold colours). This representation is similar to a heat map. It is an easy way to identify a "hot spot" (i.e. a configuration with an excessively long evacuation time) or, to the contrary, a more comfortable configuration with an acceptable evacuation duration time, which would save more lives.

In Fig. 5.18, the parameters that model the topological characteristics of the street linking node 1 to node 2, and node 2 to node 3, are set to : One can see the necessary configurations that are required for a given evacuation time in Figs. 5.18 and 5.19. As an example, if we consider a duration between 20 and 25 minutes, in Fig. 5.18, this duration can be obtained with a narrow intermediate place that can accommodate 40 to 50 people, only if the capacity of place 3 is greater than 30,000 people. Instead the same duration is possible with a smaller place 3 (with a capacity between 16,500 and 17,000 people) if the capacity of place 2 is higher (between 90 and 100 people). That means that if planners who organize sites to host public events or festivals cannot enlarge the intermediate place (for example, due to the shape of a historic city center) they must establish a larger final evacuation shelter.

Conclusion

In this article we have developed a new model of weighted human behaviour coupled on street and place networks, in the context of an urban terrorist attack, thus improving the PCR model [START_REF] Provitolo | Les comportements humains en situation de catastrophe: de l'observation à la modélisation conceptuelle et mathématique[END_REF][START_REF] Verdière | Mathematical Modeling of Human Behaviours During Catastrophic Events[END_REF] with bottleneck coupling and by taking into account the capacity of every place and the number of people stranded in these places. The simulation results in a simple network with three nodes (places or public squares) and two edges (streets) that demonstrates the key role of the capacity of an intermediate place in the dynamics of evacuation from dangerous to safe places. This model is presented with a user-friendly graphical representation, which allows planners to accurately consider where to host public events in a specific territorial context.

Fig. 5 . 1

 51 Fig. 5.1 Spatial configuration of Jallianwalla Bagh (Amritsar, India) © R. Lozi, Feb. 2018

Fig. 5 . 2

 52 Fig. 5.2 Jallianwana Bagh Memorial (Amritsar, India) © D. Provitolo, Feb. 2018 More recently, two terrorist attacks in Mediterranean cities, in the heart of the historic town centres of Barcelona (2017; Fig. 5.4) and Nice (2016), have shown that the population fled through the labyrinth of alleys to find refuge in urban places, such as public squares,

Fig. 5 . 3 Fig. 5 . 4

 5354 Fig. 5.3 Bullet marks can be easily seen on the wall (Amritsar, India) © D. Provitolo, Feb. 2018

Fig. 5 . 5

 55 Fig. 5.5 Graphical representation of the Panic, Control, Reflex behaviour model (PCR) in the exceptional situation of disaster. Source: from [3]

  .15) (Fig. 5.6).

Fig. 5 . 6 Following [ 6 ]

 566 Fig. 5.6 Forcing functions () t  , blue curve, and () t  , red curve,

Fig. 5 . 7

 57 Fig. 5.7 Possible paths of rushing people in a city network.In this article, we define the city by a mathematical graph, where the places and public squares are called vertices or nodes and are denoted 12 , , , p

Firstly

  PCR model, into the Weighted PCR (WPCR) model, by introducing new data, with Weight standing for the relative sizes of crowd, places and streets.

Fig. 5 . 8

 58 Fig. 5.8 Two nodes network: people in every state are rushing from node 1 towards node 2, maintaining the same behavioural class. (This representative colour scheme is used for the remainder of the article.)

,

  course, it is easy to relax such limitations, which are not dependent upon the coupling, by not eliminating the corresponding terms in the equations.

Fig. 5 . 9

 59 Fig. 5.9 Relationship between flow and density of pedestrians from the literature, from [17, derived from Fruin, Weidmann, Virkler, Older, Sarkar and Tanariboon]

Fig. 5 .

 5 Fig. 5.11 Graph of the bottleneck coupling function

  instead of parameters 1 B and 2 B , becomes important for the pace at which the "reservoir" of people in reflex behaviour is emptied.

Fig. 5 .

 5 Fig. 5.13 Situation 2: Convergence towards the fixed point (0, 225, 75, 0, 375, 125)

  are forbidden. People in this stunned state can only change their behaviour (from r to p or c).

Fig. 5 ..

 5 Fig. 5.15 Three node network: people in the reflex behaviour state are stranded in their original node To achieve this goal Eq. (5.29) is simply modified to vanish parameters ,, , 1, 2, 2, 3 r i j i j

Fig. 5 .

 5 Fig.5.[START_REF] Vanumu | Fundamental diagrams of pedestrian flow characteristics: A review[END_REF] Convergence towards the fixed point ( 0, 0, 0, 0, 0, 0, 0, 15000, 5000 )

Fig. 5 .

 5 Fig. 5.17 Convergence towards the fixed point ( 0, 0, 0, 0, 0, 0, 0, 15000, 5000 ) for the values of parameters = 0.2, = 0.4, = 0.3, = 0.1, ,1,2

3 W 2 :

 32 50, 100, 200, 1000 W 0). It is thus a population that is mainly in a state of reasoned behaviour that arrives in place 3 (Figs. 5.17g, h).

2 ,

 2 contrast, in Fig.5.19 their values are doubled: ,1,larger streets. In this case, shorter evacuation times are obtained.

Fig. 5 .Fig. 5 .

 55 Fig. 5.18 Time required to evacuate 80% of population from place 1 according to W2 and W3 with η = 0.0025 andV1,0 = 20 000

  Therefore, one can suppose that a terror attack is shaped by two characteristic times: s

											t
	(for start) and	e t (for end) with 0 se t t t for which
						    	( ) 1 , ( ) 0 , t t =   t t t t s e =  	(5.13)
	As an example for	I	trans	=		 2.5, 42.5	, these functions can be defined by
				0		if		0	x	42.5
		( ) t		cos 2	x	2.5 160	2	,	if	42.5	x	82.5
				1,		if x	82.5

https://geoazur.oca.eu/en/research-geoazur/2158-com2sica-how-to-comprehend-and-simulate-humanbehaviors-in-areas-facing-natural-disasters

Acknowledgements

This work has been supported by the French government, through i) the National Research Agency (ANR) under the Societal Challenge 9 "Freedom and security of Europe, its citizens and residents" with the reference number ANR-17-CE39-0008, co-financed by French Defence Procurement Agency (DGA) and The General Secretariat for Defence and National Security (SGDSN), and ii) the UCA JEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

We propose to model this bottleneck effect by a nonlinearity of a logistic type to keep the same philosophy as the authors cited above. Thus, the bottleneck coupling corresponding to Fig. 5.8, is given by the system

In this system the bottleneck coupling, concerning, for example, the controlled population that is moving from node 1 to node 2 is given by the term:

in the second equation of ( 16)

In this bottleneck coupling (5.25), ,1,2 c  is the parameter which models the topological characteristic of the street linking node 1 to node 2. The second factor 1 () ctof (5.25) reflects the pressure of controlled people in node 1 willing to escape towards node 2 and also the proportionality of people escaping with respect to people staying in node 1. Finally, the factor

shows the counter-pressure which is maximum (i.e. the term vanishes) when 2

because, in this case, there is no more room for people coming from node 1. This bottleneck coupling is nonlinear as shown on Fig. 5.11.

As we consider only transitional dynamics where ( ) In situation 2 (Fig. 5.13) when 1,0 2,0 2 700, 100, 500

, one obtains the following convergence towards the fixed point 5.13 that V1(t) tends to 300 = 800 -500. These first results lead us to consider more complicated spatial configuration networks in the following section.

Influence of the spatial configuration on the pace of evacuation

We now seek to identify the obstacles that slow the escape of the crowd in the aftermath of the initial shock by analysing the topology of the network of streets and places in a city. This information can potentially be used to improve the design of a city to facilitate the escape of a crowd towards more secure places.

For the sake of simplicity, we consider first a simplified oriented network with only three nodes

Equation of the three node network

The corresponding equation of the three nodes network is

r t r t W c t p t r t c t B r t C p t C c t c t W c t p t r t p t B r t C p t C c t p t W c t p t r t r t B B r

t W c t p t r t r t W c t p t r t c t B r t C p t C c t c t W c t p t r t c t W c t p t r t p t B r t C p t C c t

,2,3 2 3 3 3 3 (5.32)

t W t c t p t r t p W c p r r t B B r t r t W c t p t r t c t B r t C p t C c t c t W c t p t r t p t B r t C p t C