<figure><img src="image_1.png" /></figure>

HAL
open science

# Towards Parameter Identification of a Behavioral Model from a Virtual Reality Experiment

Nathalie Verdière, Oscar Navarro, Aude Naud, Alexandre Berred, Damienne

Provitolo

► To cite this version:

Nathalie Verdière, Oscar Navarro, Aude Naud, Alexandre Berred, Damienne Provitolo. Towards Parameter Identification of a Behavioral Model from a Virtual Reality Experiment. Mathematics , 2021, 10.3390/math9243175 . hal-03473102

HAL Id: hal-03473102

https://hal.science/hal-03473102

Submitted on 9 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Towards Parameter Identification of a Behavioral Model from a Virtual Reality Experiment

Nathalie Verdi ère <fnref n="1" />, Oscar Navarro <fnref n="2" />,*, Aude Naud <fnref n="2" />, Alexandre Berred <fnref n="1" /> and Damienne Provitolo <fnref n="3" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Verdière, N.; Navarro, O.; Naud, A.; Berred, A.; Provitolo, D. Towards Parameter Identification of a Behavioral Model from a Virtual Reality Experiment. Mathematics **2021**, *9*, 3175. https://doi.org/10.3390/math9243175

Academic Editor: António M. Lopes

Received: 15 November 2021

Accepted: 3 December 2021

Published: 9 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In this paper, we investigate the calibration of a mathematical model describing different behaviors occurring during a natural, a societal, or a technological catastrophe. This model was developed in collaboration with geographers and psychologists. To collect information on the level of stress, psychologists of the LPPL laboratory of Nantes (France) led virtual reality experiments. These experiments consisted in immersing individuals in a situation of catastrophe and measuring their electrocardiogram. From the physical and biological data collected, we present the methodology to calibrate the behavioral model. First, a theoretical analysis is carried out to determine (i) if the parameters can be uniquely estimated, (ii) the minimal number of discrete measurements required for the estimation. Then, from these analyses, an estimation procedure is performed to calibrate the mathematical model or at least to have an order magnitude of the model parameters. Through this work, we will show from simulations that the proposed system makes it possible to apprehend non observable human processes.

**Keywords:** parameter estimation; virtual reality; identifiability; behavioral model

## 1. Introduction

The past forty years have witnessed a significant increase in the number of natural or technological disasters. A better understanding of their consequences is essential to ensure the safety and security of populations. Improvements are particularly needed in our understanding and ability to anticipate human individual and collective behavior with respect to complex disasters of any origin [1–13]. The way we predict or anticipate human responses determines how we handle emergencies.

A few previous works concerning the modeling of human emotional behavior during a catastrophic event can be found in the literature. However, most of them focus on specific situations or phenomena, such as the feeling of panic, the fear propagation mechanism [3] or crowd behavior [14]. These works try to capture the dynamics of these mechanisms from the observation of interactions, such as on video. Panic is the emotion most studied since it is the most observable. The main difficulty, however, is to observe the different emotions of individuals during such events. Information is limited on the true evolution of a population's response during a disaster because of the difficulty of observing and analyzing human reactions in real time. The main sources of information on the behaviors adopted during a disaster are interviews and surveys carried out with operational actors, residents and victims after (or before) the disaster and in a specific region.

<footnote n="1">¹ Laboratoire de Mathématiques Appliquées, FR-CNRS-3335, 25, Rue Philippe Lebon, 76063 Le Havre, France; nathalie.verdiere@univ-lehavre.fr (N.V.); alexandre.berred@univ-lehavre.fr (A.B.)</footnote>

<footnote n="2">² Laboratoire de Psychologie des Pays de la Loire—UPRES EA 4638, University of Nantes, 44312 Nantes, France; Aude.Naud@univ-nantes.fr</footnote>

<footnote n="3">³ Damienne Provitolo, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Géoazur, Campus Azur 250, Rue Albert Einstein, CEDEX, 06905 Sophia Antipolis, France; damienne.provitolo@geoazur.unice.fr</footnote>

<footnote>* Correspondence: oscar.navarro_carrascal@unimes.fr</footnote>
Additionally, it has long been accepted by geographers and psychologists that an individual can adopt several emotional behaviors and can change behavior several times during a catastrophe leading to different group dynamics. In [15], the authors propose a model based on emotional contagion which associates the decision to search for peers with, a set of emotions. However, it was not calibrated by real experiments. It was validated from simulations showing that, by choosing certain parameters, certain types of searching decisions could be determined.

To further understand human behavior and its transitions during disasters, an inter-disciplinary study was carried out by mathematicians and researchers in the humanities and social sciences associated with field and political workers (operational and institu-tional experts) within the framework of the Com2SiCa project (ANR project Com2SiCa (https://www.com2sica.cnrs.fr/, accessed on 1 December 2021)).

The first mathematical model developed by this research group was the Reflex, Con-trol, Panic model (RCP) [5,6]. This model was based on the neurobiological response to a threat, or the stress response, with regard to the areas of the brain involved in the be-havioral responses [16]. Two categories of behaviors were proposed: instinctive behaviors commanded by the reptilian brain, and acquired and controlled behaviors commanded by the prefrontal cortex. The reflex behaviors of instinctive flight, panic, and alert are classified in the first category. The controlled behaviors are classified in the second [6]. This RCP model was then enriched by integrating knowledge from the psychology of emotions. This change in point of view has led to a new Alert-Control-Panic (APC) model integrating other emotions. It provides a better interpretation of states and behavioral state changes observed in real or simulated disaster situations. These states are not only qualified, with regard to the brain areas responding to stress, but also, with regard to the management of emotions. Two key variables were retained: the emotional load (strong or weak) and the emotional regulation (strong or weak).

To better understand the different behaviors experienced during a catastrophic event and study them in ideal conditions, a virtual reality experiment was conducted in the laboratory of Psychology of Nantes between September and November in 2020. This virtual reality experiment consisted in immersing the volunteer in a sudden, unforeseen disaster without pre-warning signs. Thanks to this new approach, observation data were obtained and allowed us to calibrate the APC model and especially to give the orders of magnitude of the parameters. As this paper will show, this calibrated model could be an interesting tool to capture and understand the different behavior dynamics that can occur during a catastrophic event.

However, two problems arose when estimating the parameters which make it difficult to implement the model in practice:

* The presence of non-linearities in the mathematical model due to the imitation processes;

* The small number of data.

This realistic but difficult context implies a prior fine analysis of feasibility of the parameter estimation, which is why we needed to address two main questions before estimating the parameters:

* Is it possible to estimate all of the model's parameters from the physiological measure-ments provided by virtual reality?

* What is the minimal number of discrete time measurements needed to be sure of the parameter estimation procedure results?

The answers were provided by two theoretical studies. The first is an identifiability analysis which ensures that, from given measurements, it is possible to uniquely estimate the model parameters [17,18]. It is an important issue when examining links between the mathematical model and data. If this step in the calibration of the model is overlooked, it can result in misleading or incorrect parameter estimates and model predictions. The sec-ond analysis determined the minimal number of time measurements required to estimate the unknown parameters.
Once these qualitative studies were conducted on the mathematical model, the procedure for estimating the parameters drawn from the data collected during the virtual reality experiment could be explained. The procedure is based on the resolution of a non-linear least-squares problem solved with the genetic algorithm [19]. A global algorithm was preferred to a local one since there was no initial guess of the parameter values.

The paper is structured as follows. Section 2 presents the different behaviors defined by psychologists and their modeling. In Section 3, the protocol implemented by the psychologists of Nantes is presented and the links between the mathematical model and the virtual reality experiment are provided in Section 4. In Section 5, the identifiability of the model is studied based on the observations made by the psychologists and the minimal number of observation points required to estimate the parameters is deduced. The parameter estimation procedure is carried out in Section 6. Section 7 concludes the paper and shows a possible application of the calibrated model to better understand the processes inducing changes in behaviors.

## 2. The Different Behaviors and Their Modeling

### 2.1. *Emotional Load and Its Regulation*

Within the framework of the Com2SiCa project, we focused on two key variables to identify and characterize different human behaviors occurring during a sudden catastrophe: the intensity of the emotional charge and its regulation. The first corresponds to the quantitative aspect of the emotion, and is adaptive to situations of threat or catastrophe [20]. We note that the intensity of the emotion is high in situations of fear or stress. The second is defined as the processes involved in adapting to relatively strong emotional episodes. Most researchers agree that emotions are not uncontrollable forces but they can be modulated when necessary (see [21] and the references within). Furthermore, adaptive emotion regulation consists in modulating—not eliminating—the experience of emotions [22,23]. By crossing these two variables, we characterized a wider range of human behaviors that may occur during a disaster than the sole panic behavior. As this behavioral diversity is difficult to model mathematically, we reduced it by identifying the behavioral states in which reactions can be located. Hence, three behavioral states were proposed:

* The alert state corresponds to a low emotional charge and low regulation. It relates to a phase of assimilation of information relative to the event, a moment when people search in a very short time for information about the scenario they are experiencing. It marks a break from everyday behavior. Alert behaviors are often described by startle movements, rapid eye movements or a fixed gaze, and verbal interrogation of people close by;

* The controlled state corresponds to an emotional charge which is more or less strong. This emotional charge is regulated which is why this behavior is controlled. People can regulate their emotions, to act and adapt their behavior to the crisis context. The controlled behavior does not ensure the safety of the person or the individual's survival. It concerns behaviors such as mutual aid, protection and care of others but also anti-social behaviors, such as theft, looting, and voyeurism;

* The panic state corresponds to a high emotional charge and a low regulation insufficient to provide a controlled behavior. Individuals in panic can adopt behaviors of panic flight, stupor, and disorderly agitation with cries of fear, for example.

Scientists in the humanities and social sciences can use different methods to differentiate these states. For example, they can associate each one with different physical manifestations that they can observe during a catastrophic event (Figure 1). They can collect physiological data to capture these behaviors and quantify them. Alert and control behaviors are associated with basic physiological states such as a low heart rate contrary to panic states. To capture these three states by simulation in a virtual reality experiment, physiological measures may be used, such as electrocardiograms and skin conductance. The aim is to quantify the emotional load and its intensity from the perspective of learned or controlled behaviors versus instinctive or automatic ones (see Section 3).
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Representation of the different behaviors and the corresponding observable physical manifestations.</figcaption></figure>

Since people do not maintain the same behaviors all the time, psychologists define two types of transitions between these three states. The first are the natural transitions specific to individuals depending on their experiences and their past. The second is imitation or social comparison which is important factor in the process of behavioral change [24]. Depending on a certain predisposition, the perception of the behavior of others encourages people to act in the same way [25]. In these collective circumstances and faced with an uncertain situation, that is to say, the absence of clear benchmarks for action, people can set up a mechanism of social comparison leading them to trust others and to consider the behavior of others as a reference and an example to follow. This form of informational imitation reduces the feeling of uncertainty.

## 2.2. Modeling the Dynamics of Human Behavior: Cross-References between the Social and Mathematical Sciences

The modeling of these different behaviors and their dynamics is presented below and is inspired by classical epidemic models, the Susceptible-Infected-Recovered (SIR) models [26]. As for SIR models, the population in situation of catastrophe is subdivided in several categories according to the adopted behavior and their different interactions are considered.

### 2.2.1. Structural Components of the APC Model

The mathematical Alert-Panic-Control (APC) model aims to describe the different behaviors during a catastrophic event, their temporal dynamics, i.e., the transitions from one state to another, and the responses to them due to qualitative parameters, i.e., the environmental context or group effects.

A catastrophic event marks a transition from everyday behavior to reactions specific to the disaster. That is why we consider a population of size *N* with an everyday behavior. We denote by *A*, *C*, and *P*, the population in the alert, the control and the panic states, respectively, and *Q* and *B* those exhibiting an everyday behavior. The first (*Q*) can be observed before the event starts but the second (*B*) occurs during or after the event. We suppose that the beginning of the catastrophe (resp. the end of the catastrophe) is modeled by the function $\gamma$ (resp. $\varphi$). Thus, when the catastrophe is triggered at time $t$ ($\gamma(t) \neq 0$), a proportion $\gamma(t)Q(t)$ of people with an everyday behavior become aware of an imminent danger and adopt an alert behavior, seeking information about what happened. Then,
depending on intrinsic and imitation processes, as explained above, people in alert state adopt either a control or a panic behavior. Finally, after a certain time, people in a control state begin to return to everyday life, by seeking shelter for example, at a proportion $\varphi(t)C(t)$ at time $t$.

Concerning the transitions (Figure 2), the parameters represented in blue, $B_i, i = 1, \dots, 4$, and $C_i, i = 1, 2$ designate the natural transitions specific to population based on their experience and past history. For $B_3$ and $B_4$, the humanities and social sciences consider that populations engaged in a panic or control behavior can return to an alert state even if the probability is low. The second type of transition, $F, G$ and $H$ in the diagram, describes the imitation or social comparison processes as explained in the following section.

<figure><img src="image_6.png" /><figcaption>Figure 2. This diagram shows the different behavior states between $Q, A, P, C, B$, and their interactions.</figcaption></figure>

### 2.2.2. Imitation Processes

Imitation processes are important transitions in threat situations. Individuals can mimic the behavior either of another individual or of individuals within a small group by observation and perception of the behavior of the other or others. In the APC model, the imitation transitions are modeled as process contamination terms in SIR models. For example, the imitation process of control behaviors by individuals in alert is given by the term $-F(A(t), C(t)) \frac{C(t)}{N} A(t)$. It models the probability for each individual with a control behavior to spread this behavior among individuals in alert and it is proportional to the portion of C in the whole population $N = Q(0) + A(0) + C(0) + P(0) + B(0)$. To express that the behavior of the majority influences the behavior of the minority through the process of imitation, we consider the following form for $F$:

$$F(A, C) = \alpha \eta \left( \frac{C}{A + \varepsilon} \right), \quad (1)$$

where $\varepsilon \ll 1, 0 \le \alpha \le 1$ and $\eta(s) = \frac{s^2}{1+s^2}$ is given in Figure 3. Since the value of the function $\eta$ depends on the proportion of people in alert and control, the probability varies according to the proportion of control and alert individuals. The function $\eta$ models that the higher the number of people with a control behavior, the higher the number of people with an alert behavior who will imitate those with a control behavior.

Analogously:

$$\left\{ \begin{aligned} G(A, P) &= \beta \eta \left( \frac{P}{A + \varepsilon} \right) \\ H(C, P) &= \gamma_{p \to c} \eta \left( \frac{C}{P + \varepsilon} \right) - \gamma_{c \to p} \eta \left( \frac{P}{C + \varepsilon} \right) \end{aligned} \right. \quad (2)$$

where $\beta, \gamma_{p \to c}, \gamma_{c \to p} \in [0; 1]$.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Representation of the function $\eta$. This sygmoidal function was chosen to model the fact that the behavior of the majority is the most imitated one.</figcaption></figure>

For the imitation term $H$, the imitation can go from panic to control and from control to panic.

Note that the imitation process can go in both directions between control and panic but not between alert/control and alert/panic, according to psychologists, since alert is a state that passes too rapidly to be imitated.

### 2.2.3. The Complete Model

The APC model comprises 5 equations, each one describing the evolution in time (in seconds in this article) of one of the population groups. Following Figure 2, the APC mathematical model is the following:

$$ 
\left\{
\begin{aligned}
\frac{dA(t)}{dt} &= \gamma(t)Q(t) - (B_1 + B_2)A(t) + B_3C(t) + B_4P(t) - F(A(t), C(t))\frac{C(t)}{N}A(t) \\
&\quad -G(A(t), P(t))\frac{P(t)}{N}A(t), \\[2ex]
\frac{dC(t)}{dt} &= B_1A(t) + C_1P(t) - (B_3 + C_2)C(t) + F(A(t), C(t))\frac{C(t)}{N}A(t) \\
&\quad +H(C(t), P(t))\frac{P(t)}{N}C(t) - \varphi(t)C(t), \\[2ex]
\frac{dP(t)}{dt} &= B_2A(t) + C_2C(t) - (B_4 + C_1)P(t) + G(A(t), P(t))\frac{P(t)}{N}A(t) \\
&\quad -H(C(t), P(t))\frac{C(t)}{N}P(t), \\[2ex]
\frac{dQ(t)}{dt} &= -\gamma(t)Q(t), \quad \frac{dB(t)}{dt} = \varphi(t)C(t).
\end{aligned}
\right.
\quad (3) $$

Model (3) is rewritten in terms of fraction of the population instead of numbers of individuals:
$$
\left\{
\begin{aligned}
\frac{da(t)}{dt} &= \gamma(t)q(t) - (B_1 + B_2)a(t) + B_3c(t) + B_4p(t) \\
&\quad - \tilde{F}(a(t), c(t))c(t)a(t) - \tilde{G}(a(t), p(t))p(t)a(t), \\[2ex]
\frac{dc(t)}{dt} &= B_1a(t) + C_1p(t) - (B_3 + C_2)c(t) + \tilde{F}(a(t), c(t))c(t)a(t) \\
&\quad + \tilde{H}(c(t), p(t))p(t)c(t) - \varphi(t)c(t), \\[2ex]
\frac{dp(t)}{dt} &= B_2a(t) + C_2c(t) - (B_4 + C_1)p(t) + \tilde{G}(a(t), p(t))p(t)a(t) \\
&\quad - \tilde{H}(c(t), p(t))c(t)p(t), \\[2ex]
\frac{dq(t)}{dt} &= -\gamma(t)q(t), \quad \frac{db(t)}{dt} = \varphi(t)c(t)
\end{aligned}
\right.
\qquad (4)
$$

with $(\tilde{\epsilon} = \frac{\epsilon}{N})$

$$
\left\{
\begin{aligned}
\tilde{F}(a, c) &= \alpha\eta\left(\frac{c}{a+\tilde{\epsilon}}\right), \quad \tilde{G}(a, p) = \beta\eta\left(\frac{p}{a+\tilde{\epsilon}}\right), \\
\tilde{H}(c, p) &= \gamma_{p\to c}\eta\left(\frac{c}{p+\tilde{\epsilon}}\right) - \gamma_{c\to p}\eta\left(\frac{p}{c+\tilde{\epsilon}}\right).
\end{aligned}
\right.
\qquad (5)
$$

It can be rewritten:

$$
\dot{x}(t, \theta) = \psi(t, x, u, \theta) \qquad (6)
$$

where $x = (a, c, p, q, b)^T = (A/N, C/N, P/N, Q/N, B/N)^T \in [0, 1]^5$, $u = (\gamma, \varphi)$, $\theta = (B_1, B_2, B_3, B_4, C_1, C_2, \alpha, \beta, \gamma_{c \to p}, \gamma_{p \to c}) \in D$ with $D := (\mathbb{R}^{+*})^6 \times (\mathbb{R}^+)^4$.

At the initial time, the entire population is supposed to be engaged in an everyday
behavior so that the following initial conditions are considered:

$$
x(0) = (0, 0, 0, 1, 0)^T. \tag{7}
$$

Following the proof of [2], the proposition is deduced:

**Proposition 1.** The Cauchy problem (6)–(7) admits a unique maximal non-negative solution for any value of the parameters $\theta \in D$.

3. The Virtual Reality (VR) Experiment

The VR experiment was conducted by psychologists in the LPPL laboratory of Nantes
(France) in a closed and controlled environment. Its aim was to determine from physio-
logical measurements, the different behavioral states and the changes occurring during a
catastrophic event by taking into account the environmental context, such as the type of
event, the spatial configuration, and the presence or absence of individuals.

3.1. *The VR Experimental Procedure and Sample*

This VR experiment consisted in immersing individuals in the simulation of a tsunami
on the Mediterranean coast of Nice (France), a scenario event which has since been included
in national regulatory tools. The scenario of the simulation was based on research data
(See [27]). To construct the spatial context, topographical information was gathered in
Nice regarding the beach, such as its width, the presence of stairs, pebbles, and access
guardrails, etc. (See [28]). Virtual characters ( avatars) were also included to observe
behavioral reactions in the context of social interactions between the survey respondents and the avatars.

Three scenarios were built based on three experimental conditions, designed to simulate different reactions, to compare them and to better understand how the context affects behavior. For example, a situation of danger does not automatically generate a feeling of threat if the person is alone. The reaction depends on the person's experience and the behaviors of those around whose behavior can be imitated. The three following scenarios were carried out:

* First scenario (neutral scenario): in this experimental condition, the person is alone without any virtual person;

* Second scenario (organized scenario): in this experimental condition, the person is in the presence of virtual persons with a suitable behavior (escape from the tsunami by walking quickly in the right direction and without shouting);

* Third scenario (non-organized scenario): in this experimental condition, the person is in the presence of virtual persons with an unsuitable behavior (escape by running in all directions and shouting).

In the second scenario, the physical reactions of the avatars express a control behavior whereas in the third one they express panic.

Eighty-eight unpaid volunteers were recruited through the local press and were assigned randomly to one of the three conditions. The sample was composed of 43 women and 45 men. Their age ranged from 18 to 56 years with a mean of 26.76 years and a standard deviation of 9.41. The participants were divided into the different conditions by respecting a balance between men and women, and a diversity in the ages. Therefore, 29 persons were assigned to the first and third scenario and 30 to the second one. The distribution of the participants with respect to their sex and the scenario is provided in Table 1.

<table><caption>**Table 1.** Distribution of the surveys with respect to their sex and the scenario.</caption><thead><tr><th rowspan="2">Condition</th><th colspan="3">Sex</th></tr><tr><th>Women</th><th>Men</th><th>Total</th></tr></thead><tbody><tr><td>First scenario</td><td>14</td><td>15</td><td>29</td></tr><tr><td>Organized scenario</td><td>15</td><td>15</td><td>30</td></tr><tr><td>Non organized scenario</td><td>14</td><td>15</td><td>29</td></tr><tr><td>Total</td><td>43</td><td>45</td><td>88</td></tr></tbody></table>

During the simulation, each participant was immersed in a virtual reality environment that provides a visual and auditory feedback (Figure 4). The subject was directly placed on the beach at the beginning of the simulation and had free movements. To capture and maintain the participant's attention, the proposed timeline was based on how a tsunami on the Mediterranean coast of Nice is perceived rather than the reality of its dynamics. As decided, after 60 s, the level of the sea decreased. Then, the wave formed on the horizon, and the water noises became louder. The wave reached the beach after 3 min. If the participants stayed immobile on the beach, the application stopped when they were struck by the wave. If the participants had an escape behavior inside the city, the application stopped when they reached a city exit. In all other cases, the application stopped 2 min after the wave submerged the beach. The steps are summed up in Figure 5.

All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of University of Angers, France (CER_20200114).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Photographs of the virtual experiment in the laboratory of Nantes. Each participant was immersed in a virtual reality environment providing visual and auditory feedback.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Process of the tsunami scenario over time: in green, progression of the wave due to the tsunami; in purple, the sounds heard by the survey respondents participating in the virtual reality experiment; in orange, the avatars action; in red, the beginning and the end of the actions.</figcaption></figure>

## 3.2. *Measurements*

Physiological measurements were recorded for each participant to determine their emotional charge at each moment of the simulation. A Nexus X MKII device was connected to the participant by electrodes. It recorded for each of them an electrocardiogram (ECG) and skin conductance in order to identify any stress signals when faced with the tsunami. The ECG is a recording that provides different variables relevant for detecting of a state of stress [29]. It reports successive polarization and depolarization of the heart, translated
into the characteristic peaks of this type of recording. In particular, it provides the Heart Rate (HR) expressed as beats per minute, calculated from the frequency of apparition of a type of peak during one minute of recording. Skin conductance measures the continuous variation in the electrical characteristics of the skin. The data are collected by applying a low, undetectable, and constant voltage to the skin and then measuring how the skin conductance varies.

In our study, only the ECG and more specifically the HR, were retained because of their reliability [30]. Even if the experiment was conducted in a laboratory, the slight movements of the participants distorted some of the conductance measurements.

Before starting each simulation, the experimenters carried out a test to check that the devices functioned properly (VR and physiological measurements) and to familiarize the participants with VR. Thus, all participants experienced the same test situation, which was an immersion in a Viking village. This test also allowed us to record a base line of the HR at rest specific to each individual. Throughout the VR experiment, the stress state of each participant was determined in comparison with the base line: a state of stress was identified as soon as the participant's HR exceeded the base line and, therefore, increased.

From the HR measurements, two categories of persons were distinguished: calm and stressful persons corresponding to a low HR and a high HR, respectively. Consequently, two series of measurements were obtained.

In the following section, we show the links between the physiological measurements and the corresponding variables of the model.

# 4. The Link between the Mathematical Model and the Experimental Data

According to [31] and Section 2, three types of behavioral states occur during a catastrophic event with respect to their emotional charge and the capacity of individuals to regulate it: the alert, control and panic states. In the following section, we associate these states with the psychological measurements.

## 4.1. *Observed Variables of the Mathematical Model*

In this section, we determine which variables or combinations (of variables) of the mathematical model can be associated with one of the experiment measures.

Based on knowledge in psychology and physiology, alert and control behaviors correspond to a basic state in which persons have a low heart rate (HR) contrary to the panic behavior (see Section 2 and [30]). During the VR experiments, two sets of measurements distinguished from the cardiac measurements were associated with a low HR or a high HR. Therefore, we consider that the first set of measurements correspond to everyday, alert and control behaviors in the APC model and the second set to the panic behavior. Persons returning to everyday life in the mathematical model can be identified by persons moving to safety in *le vieux Nice* or finishing the simulation. Their percentage during the experiment time corresponds to the difference between the total number of survey respondents and the number of control and panic respondents.

The data provided by the psychologists are gathered in Table 2 and represent the percentages of survey respondents exhibiting daily-alert-control behavior, panic behavior and return to everyday behavior in the three scenarios (neutral, organized, and non-organized). A graphical representation of the distribution of these behavioral states over time is provided in Figure 6 for the neutral scenario.

The next step was to convert the data in Table 2 into densities for the mathematical treatment. We denote $S_m$ as the density of persons with a low HR corresponding to the first set of measurements, $p_m$ the density of persons with a high HR corresponding to the second set and $b_m := 1 - S_m - p_m$ the density of persons who finished the VR scenario.

The corresponding variables of the mathematical model (4) are denoted $y_1 = q + a + c$, $y_2 = p$ and $y_3 = b$. They are also called outputs of the model. The aim of the following step was to search for the best parameters to make the trajectories of the model fit as closely as possible to the experimental data such as those of Figure 6 for the neutral scenario.
<table><caption>**Table 2.** Tabular giving the percentage of survey respondents with an everyday-alert-control behavior (EAC), panic behavior (P) and return to back to daily behavior (BD) for the three scenarios (neutral, organized, and non-organized) with respect to the time in seconds.</caption><thead><tr><th></th><th>Times (s)</th><th>40</th><th>60</th><th>80</th><th>0.100</th><th>120</th><th>140</th><th>160</th><th>180</th><th>200</th><th>220</th><th>240</th><th>260</th><th>280</th></tr></thead><tbody><tr><td rowspan="3">Neutral</td><td>P (%)</td><td>0.00</td><td>47.83</td><td>56.52</td><td>56.52</td><td>43.48</td><td>30.43</td><td>8.70</td><td>13.04</td><td>13.04</td><td>8.70</td><td>4.35</td><td>4.35</td><td>0.00</td></tr><tr><td>EAC (%)</td><td>100.00</td><td>52.17</td><td>43.48</td><td>43.48</td><td>26.09</td><td>30.43</td><td>30.43</td><td>8.70</td><td>4.35</td><td>4.35</td><td>4.35</td><td>4.35</td><td>0.00</td></tr><tr><td>BD (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>30.43</td><td>39.13</td><td>60.87</td><td>78.26</td><td>82.61</td><td>86.96</td><td>91.30</td><td>91.30</td><td>100.00</td></tr><tr><td rowspan="3">Organized</td><td>P (%)</td><td>0.00</td><td>44.00</td><td>52.00</td><td>56.00</td><td>52.00</td><td>36.00</td><td>12.00</td><td>8.00</td><td>8.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>EAC (%)</td><td>100.00</td><td>56.00</td><td>48.00</td><td>44.00</td><td>44.00</td><td>28.00</td><td>20.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>4.00</td><td>0.00</td></tr><tr><td>BD (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>4.00</td><td>36.00</td><td>68.00</td><td>88.00</td><td>88.00</td><td>96.00</td><td>96.00</td><td>96.00</td><td>100.00</td></tr><tr><td rowspan="3">Non<br>Organized</td><td>P (%)</td><td>0.00</td><td>53.85</td><td>61.54</td><td>50.00</td><td>65.38</td><td>30.77</td><td>11.54</td><td>7.69</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>EAC (%)</td><td>100.00</td><td>46.15</td><td>38.46</td><td>50.00</td><td>26.92</td><td>26.92</td><td>30.77</td><td>7.69</td><td>3.85</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>BD (%)</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>7.69</td><td>42.31</td><td>57.69</td><td>84.62</td><td>96.15</td><td>100.00</td><td>100.00</td><td>100.00</td><td>100.00</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Neutral scenario: Representation of the distribution (%) over time of the everyday-alert-control behaviors, panic behavior and return to everyday behavior.</figcaption></figure>

## 4.2. Modeling the Beginning and the End of the Catastrophe

From Figure 5 and Table 2, a simple representation of functions $\gamma$ and $\varphi$ can be found. In the scenario of the VR experiment, the survey respondents are supposed to be engaged in an everyday behavior until the 60 s point, the moment at which the beginning of the wave can be heard. Thus, $\gamma$ that models the catastrophe is supposed to be equal to 0 up to 60 s. Then, at 90 s, when the siren sounds, we suppose that all the persons become aware of the imminent arrival of danger. Consequently, 100% of survey respondents have adopted in one of the three behaviors: alert, control or panic behavior. Between 60 and 90 s, we assume that the percentage of people aware that a disaster is regularly increasing. Using the function

$$ \phi(s) = \begin{cases} 0 & \text{if } s < \tau_0 \\ 1 & \text{if } s > \tau_1 \\ \frac{1}{2} - \frac{1}{2} \cos\left(\frac{s-a}{b-a}\pi\right) & \text{if } \tau_0 \le s \le \tau_1 \end{cases} \qquad (8) $$

an approximation of $\gamma$ with $\tau_0 = 60$, $\tau_1 = 90$ is obtained (see Figure 7).

In the same way, the function $\varphi$ describing the return to an everyday behavior is constructed from $\varphi$ with $\tau_0$ (resp. $\tau_1$) as the time at which the first survey respondent (resp. the last one) left the simulation. $\varphi$ must be calculated for each scenario since persons do not adopt the same behaviors. For example, the organized scenario ends at 280 s contrary to the non-organized scenario that ends at 200 s. This difference shows that imitation processes were different depending on the scenario.
<figure><img src="image_11.png" /><figcaption>Figure 7. Representation of the function $\gamma$ that models the catastrophe using the function (8) with $\tau_0 = 60$, $\tau_1 = 90$.</figcaption></figure>

## 5. Identifiability and Minimal Number of Observation Data

This section describes the strategy implemented to estimate the parameters from the set of measurements provided in Table 2. Two strategies are possible. The first is to follow the three VR scenarios and estimate first the intrinsic parameters and then, from their knowledge, estimate the imitation processes in the organized and non-organized scenario. A second possible strategy is to estimate all the parameters from the two scenarios integrating imitation processes (organized and non-organized scenarios) without any knowledge of the intrinsic parameters.

To determine the best strategy, a global identifiability study is first carried out on the Model (6) to which the model outputs are added. The global identifiability of the model makes it possible to determine whether without noisy data and with sufficient data, all the parameters or groups of parameters can be uniquely characterized from the model outputs in all the search domain. A further specific study serves to determine whether the number of measurements is sufficient to estimate all or groups of parameters. This mathematical study will show that the second strategy does not allow us to estimate all the parameters from the second and third scenarios due to a lack of measurements as opposed to the first strategy. Therefore, the strategy selected has to correspond to VR protocol which is to first estimate the intrinsic parameters and then the imitation processes.

For the mathematical studies, we add the outputs to Equation (6) and set:

$$ \begin{cases} \dot{x}(t, \theta) = \psi(t, x(t, \theta), u(t), \theta), \\ y(t, \theta) = h(x(t, \theta), \theta) \in \mathbb{R}^3. \end{cases} \qquad (9) $$

### 5.1. Identifiability

#### 5.1.1. Definition

The identifiability property ensures that two different parameter vectors generate two distinct trajectories. More precisely:

**Definition 1.** Let $\mathcal{P}$ be an admissible set of parameters. The model (9) is said to be globally identifiable at $\theta$ if a time $t_1 \in [0, T]$ exists such that for all $t \in [0, t_1]$, for all $\bar{\theta} \in \mathcal{P}$, $\bar{\theta} \neq \theta$, the trajectories $y(t, \theta)$ and $y(t, \bar{\theta})$ are different.

#### 5.1.2. Method

To test the identifiability of the model, the method proposed in [32] and performed in [33] is also used. It is based on differential polynomials that is, expressions constructed from variables $t$, $x$, $u$, $y$, parameters $\theta_i$ and constants. Specific output-polynomials depending only on the model outputs $y$, the inputs $u$, and the parameters $\theta_i$ are obtained with the Rosenfeld-Groebner algorithm, an elimination algorithm. From the specific order
$[y, u] \prec [x]$ that consists in eliminating the unknown state variables (here the block $[v]$ means that $v$ and its derivatives are ordered by the ranking $<$ such that $v < \dot{v} < \ddot{v} < \dots$ and $\prec$ indicates a block elimination ranking in the sense "ranks lower than" between two blocks), it triangulates the differential equation system (9) using the operations of addition, subtraction, multiplication and derivation. A general set of differential polynomials with pairwise distinct leading derivatives is obtained and its solutions are the general solutions of (9) and reciprocally. The specific output-polynomials belong to this set. They can be rewritten in the following form (see [32]):

$$ P_i(y, u, \theta) = m_0^i(y, u) + \sum_{j=1}^{q_i} \Theta_j^i(\theta) m_{j,i}(y, u) = 0, \quad i \text{ from } 1 \text{ to } 3 \quad (10) $$

where $(\Theta_j^i)_{1 \le j \le q_i}$ are rational in $\theta$, $\Theta_v^i \ne \Theta_w^i$ ($v \ne w$), $(m_{j,i})_{1 \le j \le q_i}$ are differential polynomials with respect to $y$, $m_0^i \ne 0$. The highest ranking derivative which appears in $P_i$ is called the leader. For example, with the elimination order $[x_1] \prec [x_2]$, the leader of the differential polynomial $p(x_1, x_2) = x_1^2 + \dot{x}_2 + x_2 + x_1 \ddot{x}_2$ is $\ddot{x}_2$.

According to [32], the algorithm returns as many output-polynomials as outputs. In our case, we have 3 output-polynomials, and the leader of the polynomial in $P_i$ corresponds to $y_i$.

From then onwards, to simplify the notations, only one output is considered such that the index $i$ is omitted. Let us denote by $s$ the highest time derivative order of $y$ in $P$.

From the coefficients of the polynomial $P$, the initial conditions, if any, and the constraints on the parameters, we can define the real-valued function $\phi$ by

$$ \phi: \theta = (\theta_1, \dots, \theta_p) \in \mathcal{P} \mapsto (\Theta_1, \dots, \Theta_q, y(0, \theta), \dots, y^{(s-1)}(0, \theta)) \quad (11) $$

where $\mathcal{P}$ is an a priori known admissible set of parameters. Under some technical conditions omitted in this paper (for more details, the interested reader can refer to [32]), the identifiability result is the following:

**Proposition 2.** We assume that the Wronskian $W(m_1, \dots, m_q) = \det(m_k(y, u), k = 1, \dots, q)$ is not identically equal to zero (To verify this assumption, it is sufficient to verify the linear independence of the $m_k(y, u)$, $k = 1, \dots, q$. In fact, the monomials are linearly independent if there exists a time point at which the Wronskian is non-zero.) The model (9) is globally identifiable if the function $\phi$ is injective on $\mathcal{P}$.

**Remark 1.** To test the identifiability of the APC model, we used the IdentifiabilityTree algorithm developed in [33] and implemented in Maple. This algorithm provides the set of key parameters leading to the identifiability of the model and can ultimately test its global identifiability. The first step is to use the Rosenfeld-Groebner algorithm to obtain the output-polynomials. From the coefficients of these differential polynomials, the second step uses IdentifiabilityTree algorithm. This algorithm tests the identifiability of a given parameter when some parameters are known. The reason is that, in some applications, the identifiability property of the mathematical model may depend on the knowledge of some key parameters. The algorithm returns a set of parameter lists, and each set gives the identifiability of one parameter with respect to the others. In practical terms, if a parameter $\theta_1$ is not identifiable, the lists gives the parameters that has to be identified to obtain the identifiability of $\theta_1$. If the model is identifiable, only one list containing all the parameters of the model is returned.

At the end of this step, we have the following assumptions $\mathcal{H}$:

1. Functions of $\gamma$ and $\phi$ are known;

2. $y_1(t) = q(t) + a(t) + c(t)$, $y_2(t) = p(t)$, $y_3(t) = b(t)$.
### 5.1.3. Identifiability of the APC Model with No Imitation

In this section, we consider the first scenario where there is no imitation. The consequences for the mathematical model are that $F = G = H = 0$ and the parameter vector that we want to identify is $\tilde{\theta} := (B_1, B_2, B_3, B_4, C_1, C_2)$.

**Remark 2.** The identifiability of the model at $\tilde{\theta}$ when $F = G = H = 0$ means that the intrinsic parameters could be estimated uniquely from the first scenario with perfect, noisy-free data. Assuming that the profile of participants is the same in the three scenarios and, therefore, also the intrinsic parameters, the imitation parameters can be estimated from the second and third scenarios if they are identifiable.

The following proposition gives the identifiability result when $\alpha = \beta = \gamma_1 = \gamma_2 = 0$, and rewritten in an equivalent way $F = G = H = 0$.

**Proposition 1.** If there is no imitation, under the assumption $\mathcal{H}$, model (9) is identifiable with respect to $\tilde{\theta}$.

**Proof of Proposition 1.** The *IdentifiabilityTree* algorithm is applied with the elimination order $[y_1, y_2, y_3] \prec [q, a, c, p, b]$.

The output-polynomials given by the Rosenfeld-Groebner algorithm are:

$$ P_1 := (-B_1 B_2 C_2 + B_1 C_2^2 - B_2^2 B_3 - B_2^2 C_2 + B_2 B_3 C_2 + B_2 C_2^2) \dot{y}_3 \\ + B_2 C_2 \dot{y}_3 \varphi + B_2 (\dot{y}_2 + \dot{y}_2 \gamma) \varphi + (B_2^2 - B_2 C_2) \dot{y}_3 \gamma \\ + (B_1 B_2 - B_1 C_2 + B_2^2 + B_2 B_4 + B_2 C_1) \dot{y}_2 \varphi + (B_1 B_2 B_4 \\ + B_1 B_2 C_1 - B_1 B_4 C_2 - B_1 C_1 C_2 + B_2^2 C_1 - B_2 C_1 C_2) y_2 \varphi \\ + (B_2 B_4 + B_2 C_1) y_2 \gamma \varphi - B_2^2 y_1 \gamma \varphi = 0, \quad (12) $$

$$ P_2 := -B_2 (\dot{y}_3 \dot{\varphi} + \ddot{y}_3 \varphi + \dot{y}_3 \varphi^2) + (B_1 C_2 + B_2 B_3 + B_2 C_2) \dot{y}_3 \varphi \\ - B_1 \dot{y}_2 \varphi^2 + (-B_1 B_4 - B_1 C_1 - B_2 C_1) y_2 \varphi^2 = 0, $$

$$ P_3 := \dot{y}_1 + \dot{y}_2 + \dot{y}_3 = 0 $$

The third polynomial is not informative on the coefficients, thus only $P_1$ and $P_2$ will be used to study the identifiability of the model. To put these polynomials in the form (10), we let $\tilde{B}_1 := \frac{B_1}{B_2}$ ($B_2$ is supposed to be non-zero) and we divide $P_1$ and $P_2$ by $B_2$. The following real-valued function is obtained $\phi(\theta) = (B_2, C_2, B_2 - C_2, B_4 + C_1, (B_2 - C_2) \tilde{B}_1 + B_2 + B_4 + C_1, (B_2 - C_2) (\tilde{B}_1 C_2 + B_3 + C_2), ((B_4 + C_1) \tilde{B}_1 + C_1) (B_2 - C_2), \tilde{B}_1, \tilde{B}_1 C_2 + B_3 + C_2, (B_4 + C_1) \tilde{B}_1 + C_1)$.

Under the set of conditions $C(\theta) = \{\tilde{B}_1 > 0, B_2 > 0, B_3 > 0, B_4 > 0, C_1 > 0, C_2 > 0\}$, the algorithm *IdentifiabilityTree* allows us to conclude that the parameters $\tilde{B}_1, B_2, B_3, B_4, C_1, C_2$ are identifiable, and consequently $B_1$ too. $\square$

### 5.1.4. Identifiability of the Imitation Parameters

Now that the parameters of the linear part are known, an identifiability study can be carried out on the set of imitation parameters that is $\tilde{\theta} := (\alpha, \beta, \gamma_{c \to p}, \gamma_{p \to c})$. We have the following proposition.

**Proposition 3.** We assume that the intrinsic parameters are known. Under the assumption $\mathcal{H}$, the model (9) is identifiable with respect to the $\tilde{\theta}$.
**Proof.** First, according to the last equation of (4), c can be supposed as an output of the model since b and φ are known according to assumption H. Let y₄ := c. The sum of the first three equations of (4) is equal to:

$$ \dot{a} + \dot{c} + \dot{p} + \dot{q} = -\gamma q - \varphi c. \quad (13) $$

Equation (13) can also be rewritten $\dot{q} + \gamma q = \dot{y}_1 + \dot{y}_2 - \varphi y_4$. This differential equation has a unique solution with the initial condition $q(0) = 1$. Furthermore, q can also be considered too as an output of the model. Let $y_5 := q$. Finally, from (13), we deduce that $y_6 := a$ is an output of the model. Consequently, the first three equations of (9) are output-polynomials.

To prove the identifiability of the imitation parameters α and β, it is sufficient to prove that the functions $\varphi_1(t) = \frac{c(t)^3}{(a(t) + \tilde{\varepsilon})^2 + c(t)^2}$, $\varphi_2(t) = \frac{p(t)^3}{(a(t) + \tilde{\varepsilon})^2 + p(t)^2}$ are linearly independent. To test this condition, we evaluate the Wronskian

$$ W(\varphi_1, \varphi_2) = \begin{vmatrix} \varphi_1(t) & \varphi_2(t) \\ \dot{\varphi}_1(t) & \dot{\varphi}_2(t) \end{vmatrix} $$

that gives

$$ W(\varphi_1, \varphi_2) = \frac{3 p(t)^2 c(t)^2 \left(A + B - \frac{2C}{3}\right)}{Den} \quad (14) $$

with

$$ \begin{align} A &:= -((a(t) + \tilde{\varepsilon})^2 + p(t)^2) p(t) \left(\frac{c(t)^2}{3} + (a(t) + \tilde{\varepsilon})^2\right) \dot{c}(t), \nonumber \\ B &:= c(t) \left(\frac{p(t)^2}{3} + (a(t) + \tilde{\varepsilon})^2\right) \left((a(t) + \tilde{\varepsilon})^2 + c(t)^2\right) \dot{p}(t), \nonumber \\ C &:= \dot{a}(t) p(t) c(t) (c(t) - p(t)) (c(t) + p(t)) (a(t) + \tilde{\varepsilon}) \nonumber \\ Den &:= (c(t)^2 + a(t)^2 + 2a(t)\tilde{\varepsilon} + \tilde{\varepsilon}^2) (a(t)^2 + 2a(t)\tilde{\varepsilon} + p(t)^2 + \tilde{\varepsilon}^2) \neq 0. \tag{15} \end{align} $$

The Wronskian is not identically equal to zero. Therefore, α and β are identifiable. We can prove in the same manner that $\gamma_{p \to c}$ and $\gamma_{c \to p}$ are identifiable. □

**Remark 3.** To prove the identifiability of the APC model with no imitation (see Section 5.1.3) we could use the same reasoning as in the previous proof. However, we need an explicit expression of the output-polynomials to determine the minimal number of measurements in Section 5.2.

### 5.1.5. Identifiability of the Complete Model

The following proposition gives the identifiability result of the complete model.

**Proposition 4.** Under assumption H, model (9) is identifiable.

**Proof.** Using the same arguments as in the proof of Proposition 3, we find that the equations of the model are output-polynomials and that all the parameters are identifiable. □

**Remark 4.** The output-polynomials deduced from the complete model (with imitation terms) and the Rosenfeld-Groebner algorithm have a higher time derivative order than P₁ and P₂ defined at (12). Because of the non-linear terms corresponding to the imitation processes, more derivation operations are necessary in the elimination procedure to obtain the corresponding output-polynomials.

### 5.2. Minimal Number of Observation Data

In [34,35], the authors proposed a method to determine the minimal number of observation points in order to obtain the local identifiability of the system. The result is based on the generation of a function $\rho = \rho(\Theta, y, \dot{y}, \dots, y^{(k)})$, where k is a non-negative integer, such that $\rho$ has continuous partial derivatives with respect to $\Theta$. The local identifiability
results from the local inversion of $\rho$ by using the implicit function theorem in a neighborhood where the function vanishes. [34] noticed that differential polynomials obtained from differential algebra methods can generate this function $\rho$. We present below an adaptation of their result with the output-polynomial (10).

To simplify the notations in this section, we consider only one output-polynomial. The index $i$ is, therefore, omitted in (10). We denote the coefficients of $P$ by $\Theta = (\Theta_k(\theta))_{k=1,...,q}$. Note that $P_1$ and $P_2$ defined at (12) can be rewritten in this form.

**Proposition 5.** Let us assume that

1. The function

$$ \rho: \mathbb{R}^q \times \mathbb{R}^m \to \mathbb{R}^q \\ (\Theta, y) \mapsto (P(\theta, y), P(\theta, y)^{(1)}, \dots, P(\theta, y)^{(q-1)}) \qquad (16) $$

is a differentiable function with respect to the first variable;

2. The function $\phi$ defined at Equation (11) is injective;

3. The Wronskian $W(m_1, \dots, m_q)$ defined at Proposition 2 is not identically equal to zero.

In order to have the local identifiability of the model, the minimum number of measurements of $y$ is $s_y + q$, where $s_y$ is the higher order derivative of $y$ in $P$.

**Proof.** First, note that by definition $\rho$ vanishes for all $y$ and $\Theta$. Then, let us apply the Taylor expansion of $\rho$ at $\Theta^*$

$$ \rho(\Theta, y) \approx \rho(\Theta^*, y) + (\Theta - \Theta^*) \frac{\partial \rho}{\partial \Theta}(\Theta^*, y). \qquad (17) $$

By definition, $\rho(\Theta^*, y) = 0$ and $\rho$ has continuous partial derivatives with respect to $\Theta$, thus

$$ \left| \frac{\partial \rho}{\partial \Theta} \right| = \begin{vmatrix} \frac{\partial P}{\partial \Theta_1} & \cdots & \frac{\partial P}{\partial \Theta_q} \\ \vdots & \ddots & \vdots \\ \frac{\partial P^{(q-1)}}{\partial \Theta_1} & \cdots & \frac{\partial P^{(q-1)}}{\partial \Theta_q} \\ m_1 & \cdots & m_q \\ \vdots & \ddots & \vdots \\ m_1^{(q-1)} & \cdots & m_q^{(q-1)} \end{vmatrix} \qquad (18) $$

By assumption, $|\frac{\partial \rho}{\partial \Theta}|$ is the Wronskian $W(m_1, \dots, m_q)$ which is supposed to not be identically equal to zero. According to (17), $\Theta$ can be expressed as a differentiable function of $y$ and its derivatives up to its $s_y + q - 1$ order derivative. Or to evaluate a derivative of order $s_y + q - 1$, $s_y + q$ points are needed. Consequently, to have the local identifiability of the model, the number of measurements of $y$ must be at least $s_y + q$. $\square$

**Remark 5.** In the case where Assumption 2. is difficult to check, i.e., the functional determinant $|\frac{\partial \rho}{\partial \Theta}|$ is difficult to calculate because it may require high formal calculus, it can be replaced by the following:

$q$ distinct time points $t_1, \dots, t_q$ exist such that the determinant

$$ \left| \frac{\partial P}{\partial \Theta_1}(t_1) \quad \dots \quad \frac{\partial P}{\partial \Theta_q}(t_1) \atop \frac{\partial P}{\partial \Theta_1}(t_q) \quad \dots \quad \frac{\partial P}{\partial \Theta_q}(t_q) \right| \qquad (19) $$

is not identically equal to zero.
**Example 1.** Consider the output-polynomials $P_1$ and $P_2$ given in (12). They are needed to identify all the parameters. After dividing $P_1$ by $B_2$, we obtain a polynomial of the form (10) composed of 7 coefficients and whose indeterminates are $y_1, y_2, y_3$. Its leader is $\ddot{y}_2$. This polynomial must be derived six times to get a square system. The higher order derivative in $P_1$ is 2 for $y_2$ and 1 for $y_3$. Consequently, the higher order derivatives of the system composed of 7 equations are 8 for $y_2$ and 7 for $y_1$. To evaluate a derivative of order 8 (resp. order 7), 9 points are needed (resp. 8 points). Finally, 9 measurements are needed to estimate the coefficients of the output-polynomial $P_1/B_2$.

Using the same argument with $P_2/B_2$, we deduce that $3 + 2 = 5$ measurements of $y_3$ (resp. $3 + 1 = 4$ measurements of $y_2$) are needed to estimate its coefficients.

In conclusion, the minimal number of measurements required to have the local identifiability is 9.

**Remark 6.**

* From the results of Section 3, we can conclude that a parameter estimation procedure can be implemented to estimate the intrinsic parameters using the data provided in Table 1. There is a unique parameter vector that allows the model outputs to fit the data of this table in the case of noisy-free measurements;

* One could be tempted to estimate all the parameters of the model from the second and third scenario. However, from Remark 4, order derivatives of output-polynomials are higher and there are not enough measurement points to have the local identifiability.

In summary, the parameter identification will be carried out in three steps corresponding to the three conditions in the VR experiment. The main assumption is that the participants have the same profile in the three scenarios. For each one, a set of parameters will be estimated with respect to each scenario. The intrinsic parameters will be first estimated for the Neutral scenario in which there is no imitation process. Then, the parameters describing control imitation processes will be calculated from the organized scenario and the parameters describing the panic imitation processes from the non-organized scenario. Figure 8 summarizes the steps of identification.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Diagram summarizing the steps of the identification of the model parameters. The natural transitions will be first estimated from the first scenario. Then, imitation processes will be identified from the second and third scenarios.</figcaption></figure>
# 6. Parameter Estimation

The parameter estimation procedure is accomplished following the steps defined in the previous section.

## 6.1. Minimization Problem

Given discrete-time measurements and a parametrized model describing a continuous-time process, we seek the parameter values best fitting the experimental data. To fulfill this objective, an optimization procedure is implemented to adjust the model parameters and the measurements. It is constructed from a $L^2$ cost function given by the sum of the square of the differences between the outputs of the model and the experimental data at each discrete time.

Let $(t_i)_{i=1,...,M}$ be the discrete time instants corresponding to the physiological measurements, $S_m, p_m$, and $b_m$, defined in Section 4.1. These three vectors correspond to the densities of persons engaged in a control, panic and return to everyday behavior at each discrete time. Set $y^{mes} = (S_m, p_m, b_m)^T$ and $y(\theta) = (y_1(t_i), y_2(t_i), y_3(t_i))_{i=1,...,M}^T$ obtained from the outputs of the model evaluated at each $t_i$ with the parameter vector $\theta$, where the exponent “T” indicates the transpose of the matrix.

The cost function $J(\theta)$ is the square of the sum of the errors, that is

$$J(\theta) = \|y^{mes} - y(\theta)\|_2^2 = \sum_{i=1}^{M} \sum_{j=1}^{3} (y_{i,j}^{mes} - y(\theta)_{i,j})^2 \quad (20)$$

where $y_{i,j}^{mes}$ (resp. $y(\theta)_{i,j}$) is the component at the i-th line and j-th column of the matrix $y^{mes}$ (resp $y(\theta)$) and M is the number of discrete-time measurements.

The solution of the optimization algorithm is denoted

$$\hat{\theta} = \operatorname*{argmin}_{\theta \in \mathcal{P}} J(\theta) \quad (21)$$

Several local and global optimization algorithms are proposed in the literature to solve this non-linear least-squares fitting [36,37]. The results of the optimization algorithms are the sets of parameters that produce the best fit between the simulations and experimental data.

In our case, local optimization algorithms such as Levenberg-Marquardt are not reliable for finding the global minimum since they require a first initial guess of the parameters that we do not have.

Global optimization methods and particularly stochastic ones such as simulated annealing, genetic algorithms, and evolutionary computation are suggested as the methods of choice when knowledge of the parameter values is lacking. Starting with a suitable initial domain, global optimization methods search more or less exhaustively through the parameter values in the attempt to minimize the cost function. One of the main problems associated with these methods is that they tend to be computationally expensive and may not perform well when the noise in the measurements is significant.

In this work, we chose the genetic algorithm inspired by Darwin's theory of biological evolution [19,38]. This evolutionary algorithm imitates genetic processes. According to predefined operations, it generates children by selecting and combining/muting the best parents in the current generation. It reiterates this process until the population evolves into an optimal solution. Since it uses probabilistic transition rules instead of deterministic ones, it can search for solutions in a large domain. Furthermore, contrary to other algorithms, function derivatives are not required. For more details, the readers can refer to [39].

In the following section, we present the numerical results.
## 6.2. Results

The function optim_ga implemented in the software Scilab was used to estimate the parameters. The model outputs of the mathematical model are not very sensitive to the imitation parameters. We therefore just want to have an order of magnitude of these parameters.

To estimate the numerical errors in the parameters, the sum of the relative errors between the measurement curves and the model outputs with respect to time is drawn for each case.

### 6.2.1. Scenario with No Avatars

In this scenario, 29 persons were surveyed but only 23 were kept because of the quality of their measurements. Indeed, either their signal was absent or it had too many artifacts due to little movements distorting the recorded measurements. In the optimization problem (21), the variables and their constraints are $\theta = (B_1, B_2, B_3, B_4, C_1, C_2) \in \mathcal{P}$ where

$$ \mathcal{P} = \{\theta \in [0.01, 0.7]^2 \times [0.01, 0.2]^2 \times [0.01, 0.7]^2 / B_2 + C_2 < 1\}. $$

Condition $B_2 + C_2 < 1$ was added to limit the panic behavior and make the algorithm convergent considering that the model is very sensitive to the intrinsic parameters and large values make the Runge-Kutta algorithm that solves system (4) divergent.

The genetic algorithm returns the following values for the intrinsic parameters: $B_1 = 0.1915$, $B_2 = 0.4374$, $B_3 = 0.1009$, $B_4 = 0.1741$, $C_1 = 0.1448$, $C_2 = 0.4676$.

The left subplot of Figure 9 represents the trajectories of the real data $S_m, p_m, b_m$, and the trajectories of the outputs calculated with the estimated parameters. In the right subplot, the curves represent the simulations of $a, p, c, b, q$.

<figure><img src="image_13.png" /><figcaption>**Figure 9. (Left)**: Representation of real data $S_m, p_m, b_m$, and the output trajectories: densities of the control ($y_1$), panic ($y_2$), and everyday ($y_3$) behavior after the catastrophe. **(Right)**: Simulations of $a, p, c, b, q$ with the parameters obtained from the genetic algorithm.</figcaption></figure>

From the three output trajectories, the sum of the relative errors between the measurement curves is calculated and plotted with respect to time in Figure 10. The relative error increases between 120 and 150 s when individuals begin to move to safety in *le vieux Nice* but slightly decreases until the end of the scenario. Consequently, the mathematical model fits well with the beginning of the measurement curves. There is a slight degradation in the results when individuals begin to stop the VR but the relative error stays constant and decreases a little when considering a longer time.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Representation of the relative error with respect to time. It corresponds to the sum of the relative errors between $S_m, p_m, b_m$ and $y_1, y_2, y_3$, respectively. We note a deterioration of the results after 120 s but the relative error stabilizes and even decreases when considering a longer time.</figcaption></figure>

**Remark 7.** We note that the population has a low risk culture. The relations $B_2 > B_1$ and $C_2 > C_1$ imply that individuals tend to engage in a panic behavior on their own.

### 6.2.2. Organized Scenario

In the second scenario, 25 persons were selected again among 29. The values of the intrinsic parameters found in the first scenario were used. In this scenario, since there are no panic avatars, we assume that $\beta = \gamma_{c \to p} = 0$.

The optimization problem is solved with $\theta = (\alpha, \gamma_{p \to c})$ and $P = [0, 9]^2$. The results obtained with the genetic algorithm are $\hat{\theta} = (\alpha, \gamma_{p \to c}) = (0.8338, 0)$. Figure 11 gives the representation of real data and the output trajectories evaluated with $\hat{\theta}$ and the simulations of the state variables of the mathematical model.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Organized scenario. (Left): Representation of real data $S_m, p_m, b_m$, and the output trajectories: densities of control ($y_1$), panic ($y_2$), and everyday ($y_3$) behavior after the catastrophe. (Right): Simulations of $a, p, c, b, q$ with the parameters obtained from the genetic algorithm.</figcaption></figure>

From the three output trajectories, the sum of the relative errors with the measurement curves is calculated and plotted with respect to time in Figure 12.

We note a peak in the relative error corresponding to the moment that people begin to exit from VR. However, the relative error decreases when considering a longer time.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Representation of the relative error with respect to time. It corresponds to the sum of the relative errors between $S_m, p_m, b_m$ and $y_1, y_2, y_3$, respectively. At 140 s, we notice a peak in the relative error corresponding to the moment that people begin to exit from VR. However, the relative error decreases when considering a longer time.</figcaption></figure>

### 6.2.3. Non-Organized Scenario

For the third scenario, 26 persons were kept for the mathematical study. As in the previous scenario, the intrinsic parameter values found in the first case were used and since there were no control avatars, we suppose that $\alpha = \gamma_{p \to c} = 0$. In the optimization problem (21), $\theta = (\beta, \gamma_{c \to p})$ and $\mathcal{P} = \{p \in [0, 1]^2 | \beta + \gamma_{c \to p} < 0.8\}$. We need to constrain the parameter values to limit the imitation processes otherwise the panic behavior is overestimated. The result obtained with the genetic algorithm is $\hat{\theta} = (0.6599, 0.0273)$. According to Figure 13, the dynamics of the different behaviors were found once again.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Non-organized scenario. (**Left**): Representation of real data $S_m, p_m, b_m$, and the output trajectories: densities of control ($y_1$), panic ($y_2$), and everyday ($y_3$) behavior after the catastrophe. (**Right**): Simulations of $a, p, c, b, q$ with the parameters obtained from the genetic algorithm.</figcaption></figure>

The sum of the relative errors with respect to time is plotted in Figure 14. In that case, the relative error increases at 140 s when people begin to finish the RV experiment and decreases when considering a longer time.
<figure><img src="image_18.png" /><figcaption>Figure 14. Representation of the relative error with respect to time. It corresponds to the sum of the relative errors between $S_m, p_m, b_m$, and $y_1, y_2, y_3$, respectively. The relative error increases at 140 s when people begin to finish the VR experiment and decreases when considering a longer time.</figcaption></figure>

## 7. Conclusions and Research Perspectives

From discussions with geographers and psychologists, a mathematical model was created which describes the different behaviors occurring during a catastrophe. It is based on two key variables used to identify and characterize the different human behaviors in such events: the intensity of the emotional charge and its regulation. This first modeling work describes the behavioral dynamics in the context of a sudden catastrophe. However, due to difficulties in recovering data, such models are rarely calibrated. In this paper, a new step in the modeling process was achieved by calibrating the mathematical model derived from a virtual reality experiment conducted by humanities and social scientists. These scientists contributed by bringing answers to non observable processes. For example, from the mathematical model, a finer study of the transition process in behavioral changes could be carried out leading to a better understanding of these processes and consequently an improvement in decision-making strategies. In the experimental conditions described in this paper, we can characterize what was the most influential process, as explained below.

To simplify the notations, we let $\tilde{H}(c, p) = \tilde{H}_{p \to c}(c, p) - \tilde{H}_{c \to p}(c, p)$ with $\tilde{H}_{p \to c}(c, p) = \gamma_{p \to c} \eta \left( \frac{c}{p+\tilde{\epsilon}} \right)$ and $\tilde{H}_{c \to p}(c, p) = \gamma_{c \to p} \eta \left( \frac{p}{c+\tilde{\epsilon}} \right)$.

We introduce the following functions given the inflow due to imitation and intrinsic processes towards control and panic behaviors with respect to the time. For $t \in [0, T]$, the process inducing control behaviors by imitation is

$$Im_c(t) = \tilde{H}_{p \to c}(c(t), p(t)) c(t) p(t) + \tilde{F}(a(t), c(t)) c(t) a(t),$$

and for the intrinsic process

$$In_c(t) = B_1 a(t) + C_1 p(t).$$

In the same way, the imitation and intrinsic processes inducing panic behaviors are, respectively,

$$Im_p(t) = \tilde{H}_{c \to p}(c(t), p(t)) c(t) p(t) + \tilde{G}(a(t), p(t)) p(t) a(t),$$

and for the intrinsic process

$$In_p(t) = B_2 a(t) + C_2 c(t).$$
The results of the organized scenario are provided in Figure 15 and of the non-organized scenario in Figure 16. We first observe that intrinsic processes seem to dominate imitation processes in both scenarios. Second, we notice that imitation processes are more important in the non-organized scenario which shows that a panic behavior is more frequently imitated. This preliminary result may confirm the benefits of training populations in risk management policies. To confirm this result, further experiments with specific assumptions should be conducted in the future.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Inflow due to imitation and intrinsic processes towards control with respect to time in the organized scenario.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Inflow due to imitation and intrinsic processes towards panic with respect to time in the non-organized scenario.</figcaption></figure>

**Author Contributions:** O.N. and A.N. designed and performed the virtual reality experiment. N.V. analyzed data, performed the mathematical analysis, and wrote the manuscript. A.B., D.P. and O.N. reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the French government, through 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).

**Institutional Review Board Statement:** The protocol has been approved by the ethics committee of the University of Angers, France (CER_20200114).
**Informed Consent Statement:** The participants were informed of the aims of the study, and their informed consent to participate was obtained. The participants were informed that they could stop the experiment whenever they wanted without giving an explanation.

**Data Availability Statement:** The datasets supporting this article are available at https://lmah.univ-lehavre.fr/~verdiere/Table_DataF, accessed on 1 October 2021. The program to test the relative identifiability is available at https://lmah.univ-lehavre.fr/~verdiere/ModeleAPC, accessed on 1 October 2021.

**Conflicts of Interest:** The authors declare no competing interests.

## References

1. Boschetti, L.; Provitolo, D.; Tric, E. A method to analyze territory rezilience to natural hazards, the example of the French Riviera against tsunami. In Proceedings of the EGU General Assembly Conference Abstracts, Vienna, Austria, 23–28 April 2017; Volume 19, p. 12935.

2. Cantin, G.; Verdière, N.; Lanza, V.; Aziz-Alaoui, M.; Charrier, R.; Bertelle, C.; Provitolo, D.; Dubos-Paillard, E. Mathematical modeling of human behaviors during catastrophic events: Stability and bifurcations. *Int. J. Bifurc. Chaos* **2016**, *26*, 1630025. [CrossRef]

3. Cornes, F.; Frank, G.; Dorso, C. Fear propagation and the evacuation dynamics. *Simul. Model. Pract. Theory* **2019**, *95*, 112–133. [CrossRef]

4. Helbing, D.; Johansson, A. Pedestrian, Crowd and Evacuation Dynamics. *Encycl. Complex. Syst. Sci.* **2010**, *16*, 697–716. [CrossRef]

5. Verdière, N.; Lanza, V.; Charrier, R.; Provitolo, D.; Dubos-Paillard, E.; Bertelle, C.; Aziz-Alaoui, M. Mathematical modeling of human behaviors during catastrophic events. In Proceedings of the International Conference on Complex Systems and Applications, Le Havre, France, 23–26 June 2014; pp. 67–74.

6. Provitolo, D.; Dubos-Paillard, E.; Muller, J.P. Emergent human behaviour during a disaster: Thematic versus complex systems approaches. In Proceedings of the European Conference on Complex System, Vienna, Austria, 15 September 2011; pp. 1–11.

7. Reghezza-Zitt, M.; Rufat, S. *Resilience Imperative: Uncertainty, Risks and Disasters*; Elsevier: Amsterdam, The Netherlands, 2015.

8. Venel, J. Mathematical and Numerical Modelling of Crowd Motion. Ph.D. Thesis, Université Paris Sud—Paris XI, Orsay, France, 2008.

9. Verdière, N.; Cantin, G.; Provitolo, D.; Lanza, V.; Dubos-Paillard, E.; Charrier, R.; Aziz-Alaoui, M.; Bertelle, C. Understanding and simulation of human behaviors in areas affected by disasters: From the observation to the conception of a mathematical model. *Glob. J. Hum. Soc. Sci.* **2015**, *15*, 7–15.

10. Lanza, V.; Charrier, R.; Verdière, N.; Dubos-Paillard, E.; Navarro, O.; Provitolo, D.; Bertelle, C.; Cantin, G.; Aziz-Alaoui, M. Mathematical and geographical approach in the modeling of a network of human behavioral systems. In Proceedings of the XTerM2019, Le Havre, France, 26–28 June 2019.

11. Wheatley, T.; Wegner, D.M. Psychology of Automaticity of Action. In *International Encyclopedia of the Social and Behavioral Sciences*; Elsevier: Amsterdam, The Netherlands, 2001.

12. Zachariadis, V.; Amos, J.; Kohn, B. *Simulating Pedestrian Route Choice Behaviour under Transient Traffic Conditions*; Emerald Group Publishing Limited: Bingley, UK, 2009; pp. 113–135. [CrossRef]

13. Zeng, Z.; Nakamura, H.; Chen, P. The 9th International Conference on Traffic and Transportation Studies (ICTTS 2014). A Modified Social Force Model for Pedestrian Behavior Simulation at Signalized Crosswalks. *Procedia-Soc. Behav. Sci.* **2014**, *138*, 521–530.

14. Mehdi, M.; Kapadia, M.; Thrash, T.; Sumner, R.; Gross, M.; Helbing, D.; Holscher, C. Crowd behaviour during high-stress evacuations in an immersive virtual environment. *J. R. Soc. Interface* **2016**, *13*, 20160414. [CrossRef]

15. Mao, Y.; Fan, Z.; Zhao, J.; Zhang, Q.; He, W. An emotional contagion based simulation for emergency evacuation peer behavior decision. *Simul. Model. Pract. Theory* **2019**, *96*, 101936. [CrossRef]

16. Selye, H. *Stress in Health and Disease*; Butterworths: Boston, MA, USA, 1976.

17. Daly, A.C.; Gavaghan, D.; Cooper, J.; Tavener, S. Inference-based assessment of parameter identifiability in nonlinear biological models. *J. R. Soc. Interface* **2018**, *15*, 20180318. [CrossRef] [PubMed]

18. Villaverde, A.F.; Tsiantis, N.; Banga, J.R. Full observability and estimation of unknown inputs, states and parameters of nonlinear biological models. *J. R. Soc. Interface* **2019**, *16*, 20190043. [CrossRef]

19. Haupt, R.L.; Haupt, S.E. *Practical Genetic Algorithms*, 2nd ed.; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2003.

20. Russell, J.A. A circumplex model of affect. *J. Personal. Soc. Psychol.* **1980**, *39*, 1161–1178. [CrossRef]

21. Gross, J.J. Emotion Regulation: Current Status and Future Prospects. *Psychol. Inq.* **2015**, *26*, 1–26. [CrossRef]

22. Gratz, K.; Roemer, L. Multidimensional Assessment of Emotion Regulation and Dysregulation: Development, Factor Structure, and Initial Validation of the Difficulties in Emotion Regulation Scale. *J. Psychopathol. Behav. Assess.* **2004**, *26*, 41–54. [CrossRef]

23. Kaufman, E.A.; Xia, M.; Fosco, G.; Yaptangco, M.; Skidmore, C.R.; Crowell, S.E. The difficulties in emotion regulation scale short form (DERS-SF): Validation and replication in adolescent and adult samples. *J. Psychopathol. Behav. Assess.* **2015**, *38*, 443–455. [CrossRef]

24. Festinger, L. A Theory of Social Comparison Processes. *Hum. Relations* **1954**, *7*, 117–140. [CrossRef]
25. Bargh, J.A.; Chen, M.; Burrows, L. Automaticity of social behavior: Direct effects of trait construct and stereotype activation on action. *J. Personal. Soc. Psychol.* **1996**, *71*, 230–244. [CrossRef]

26. Gao, L.Q.; Hethcote, H.W. Disease transmission models with density-dependent demographics. *J. Math. Biol.* **1992**, *30*, 717–731. [CrossRef] [PubMed]

27. Boschetti, L.; Mansour, I.; Provitolo, D.; Tric, E.; Grilli, S.; Nemati, F.; Larroque, C. The Mediterranean French coastal exposed to tsunami risk: From hazard to territorial vulnerability studies. In Proceedings of the AAG 2019, Washington, DC, USA, 3–7 April 2019.

28. Verhulst, E.; Richard, P.; Provitolo, D.; Navarro, O. Virtual Tsunami: Navigation technique and user behavior analysis during emergency. In Proceedings of the 1st International Conference for Multi-Area Simulation—ICMASim, Angers, France, 8–10 October 2019.

29. Laborde, S.; Mosley, E.; Thayer, J. Heart Rate Variability and Cardiac Vagal Tone in Psychophysiological Research—Recommendations for Experiment Planning, Data Analysis, and Data Reporting. *Front. Psychol.* **2017**, *8*, 213. [CrossRef] [PubMed]

30. Boonnithi, S.; Phongsuphap, S. Comparison of heart rate variability measures for mental stress detection. In Proceedings of the 2011 Computing in Cardiology, Hangzhou, China, 18–21 September 2011; Volume 38.

31. Navarro, O.; Krien, N.; Rommel, D.; Deledalle, A.; Lemée, C.; Coquet, M.; Mercier, D.; Fleury-Bahi, G. Coping strategies regarding coastal flooding risk in a context of climate change in a French Caribbean island. *Environ. Behav.* **2020**, *53*, 636–660. [CrossRef]

32. Denis-Vidal, L.; Joly-Blanchard, G.; Noiret, C.; Petitot, M. An algorithm to test identifiability of non-linear systems. In Proceedings of the 5th IFAC NOLCOS, St. Petersburg, Russia, 4–6 July 2001; Volume 7, pp. 174–178.

33. Verdière, N.; Orange, S. A systematic approach for doing an a priori identifiability study of dynamical nonlinear models. *Math. Biosci.* **2019**, *308*, 105–113. [CrossRef]

34. Miao, H.; Xia, X.; Perelson, A.S.; Wu, H. On identifiability of nonlinear ODE models and applications in viral dynamics. *SIAM Rev.* **2011**, *53*, 3–39. [CrossRef] [PubMed]

35. Xia, X.; Moog, C.H. Identifiability of nonlinear systems with application to HIV/AIDS models. *IEEE Trans. Autom. Control* **2003**, *48*, 330–336. [CrossRef]

36. Kiranyaz, S.; Ince, T.; Gabbouj, M. *Optimization Techniques: An Overview*; Springer: Berlin/Heidelberg, Germany, 2014.

37. Walter, E.; Pronzato, L. *Identification of Parametric Models from Experimental Data*; Springer: Berlin/Heidelberg, Germany, 1997.

38. Kazarlis, S.A.; Papadakis, S.E.; Theocharis, J.B.; Petridis, V. Microgenetic algorithms as generalized hill-climbing operators for GA optimization. *IEEE Trans. Evol. Comput.* **2001**, *5*, 204–217. [CrossRef]

39. Gen, M.; Cheng, R. *Genetic Algorithms and Engineering Optimization*; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1999.