*2.1. PPP Positioning Based on PPP-B2b Service*

As is known, precise satellite orbit and clock products play a key role in PPP processing. For PPP-B2b service, the satellite orbit corrections and clock offset corrections are broadcast by BDS-3 GEO satellites. Up to now, the PPP-B2b service has provided corrections only for BDS-3/GPS satellites. According to the PPP-B2b Interface Control Document (ICD) from China Satellite Navigation Office (CSNO) [26], the PPP-B2b orbit corrections are given in the radial (*δOradial*), along-track - *δOalong* , and cross-track (*δOcross*) directions. Hence, the orbit correction vector *δOB*2*b* = *δOradial δOalong δOcrossT* should first be transformed to the Earth-Center Earth-Fixed (ECEF) frame. This is because the satellite positions derived from CNAV1/LNAV broadcast ephemeris are based on the ECEF frame. This transformation can be described as follows:

$$
\delta \mathbf{X}^{sat} = \begin{bmatrix} \mathbf{e}\_{radial} & \mathbf{e}\_{along} & \mathbf{e}\_{cross} \end{bmatrix} \cdot \delta \mathbf{O}\_{B2b} \tag{1}
$$

with

$$\begin{array}{c} \mathbf{e}\_{radial} = \frac{r}{|r|}\\ \mathbf{e}\_{cross} = \frac{r \times \bar{r}}{|r \times \bar{r}|}\\ \mathbf{e}\_{along} = \mathbf{e}\_{cross} \times \mathbf{e}\_{radial} \end{array} \tag{2}$$

where *δXsat* = *δOx δOy δOz T* represents the PPP-B2b orbit correction vector in the ECEF frame. *r* and . *r* are the satellite position and velocity vectors derived from broadcast ephemeris. Then, the precise satellite position vector *Xprec* can be computed by

⎧ ⎪⎪⎨ ⎪⎪⎩

$$\mathbf{X}\_{\text{precc}} = \mathbf{X}\_{\text{brdc}} - \delta \mathbf{X}^{\text{sat}} \tag{3}$$

where *Xbrdc* is the satellite position vector derived from broadcast ephemeris. As described in ICD document, the PPP-B2b clock correction parameter is defined as offset to the broadcast ephemeris clock in meters. The precise satellite clock offset is given by

$$dt\_{proc}^{sat} = dt\_{brdc}^{sat} - \frac{\mathbf{C}\_0}{c} \tag{4}$$

where *C*0 represents the PPP-B2b clock correction parameter; *dtsat brdc* is the satellite clock offset derived from broadcast ephemeris; *dtsat prec* denotes the precise PPP-B2b clock offset; *c* is the velocity of light in a vacuum.

The GNSS raw code and carrier-phase measurements can read as [27,28]

$$\begin{cases} \begin{aligned} P\_i &= \rho + c \cdot (dt\_r - dt^s) + T + \frac{f\_1^2}{f\_i^2} \cdot I\_1 + B\_{r,i} - B\_i^s + \varepsilon\_{P\_i} \\ L\_i &= \rho + c \cdot (dt\_r - dt^s) + T - \frac{f\_1^2}{f\_i^2} \cdot I\_1 + \lambda\_i \left( N\_i + b\_{r,i} - b\_i^s \right) + \varepsilon\_{L\_i} \end{aligned} \end{cases} \tag{5}$$

where *i* represents the frequency number; *Pi* and *Li* are raw code and phase measurements; *ρ* denotes the geometric distance from satellite to receiver; *dtr* and *dts* are the clock offsets at the receiver and satellite, respectively; *fi* is the frequency value, *T* is the tropospheric delay, and *I*1 denotes the ionospheric delay for *L*1; *Br*,*i* and *Bs i* represent the code bias at the receiver/satellite end; *br*,*i* and *bs i* denote the phase bias at the receiver/satellite end; *εPi* and *εLi* are the unmodelled errors of code/phase measurements. The relativistic effect, Sagnac effect, Shapiro time delay [29], site displacements [30], and phase windup [31] should be corrected according to the corresponding models.

The ionosphere-free code and phase combinations are usually adopted by PPP to eliminate the first-order ionospheric delay. The ionospheric-free code and phase combinations read as [32]:

$$\begin{cases} \begin{aligned} P\_{IF} &= \mathfrak{a} \cdot P\_{\tilde{t}} + (1 - \mathfrak{a}) \cdot P\_{\tilde{j}} = \mathfrak{\rho} + \mathfrak{c} \cdot d\tilde{t}\_{\tilde{r}} - \mathfrak{c} \cdot dt^{s} - B\_{IF}^{s} + T + \mathfrak{c}\_{P\_{IF}} \\ L\_{IF} &= \mathfrak{a} \cdot L\_{\tilde{i}} + (1 - \mathfrak{a}) \cdot L\_{\tilde{j}} = \mathfrak{\rho} + \mathfrak{c} \cdot d\tilde{t}\_{\tilde{r}} - \mathfrak{c} \cdot dt^{s} + T + \lambda\_{IF} \cdot \overline{\mathcal{N}}\_{IF} + \mathfrak{c}\_{L\_{IF}} \end{aligned} \end{cases} \tag{6}$$

where *α* = *f* 2 *i* / - *f* 2 *i* − *f* 2 *j* and *λIF*·*NIF* = *λIFNIF* − *bs IF* + *br*,*IF* − *Br*,*IF*, *bs IF* and *br*,*IF* are phase bias of ionospheric-free phase combination at the satellite/receiver end, *Bs IF* and *Br*,*IF* represent code bias of ionospheric-free code combination at the satellite/receiver end, *λIF* and *NIF* are the ionospheric-free wavelength and the ionospheric-free ambiguity, respectively; *dtr* is the recombined receiver clock offset, which absorbs the receiver code bias of ionospheric-free code combination; *εPIF* and *εLIF* denote the unmodelled errors of

ionospheric-free code and phase combinations. It is noted that the ionospheric-free code bias at the satellite end should be corrected by applying PPP-B2b differential code bias corrections [15].

When the recovered PPP-B2b precise satellite orbits and clock offsets have been applied, the satellite orbit and clock errors are considered eliminated. The tropospheric delay can be divided into the dry and wet parts. The Saastamoinen model is usually used to correct the dry part, and the wet part must be estimated as unknown. Then the ionospheric-free code and phase combinations can be linearized as

$$\begin{cases} \begin{aligned} p\_{IF} &= -\mathbf{e} \cdot \mathbf{x} + \mathbf{c} \cdot d\overline{t}\_{I} + M\_{W} \cdot zwd + \varepsilon\_{P\_{IF}}\\ l\_{IF} &= -\mathbf{e} \cdot \mathbf{x} + \mathbf{c} \cdot d\overline{t}\_{I} + M\_{W} \cdot zwd + \lambda\_{IF} \cdot \overline{N}\_{IF} + \varepsilon\_{L\_{IF}} \end{aligned} \end{cases} \tag{7}$$

where *e*, *x* represent the unit vector from receiver to satellite, and the vector of position increment, respectively; *zwd* denotes the zenith wet delay and *MW* is the corresponding mapping function. In this equation, the remaining unknown parameters include only the receiver position increment vector *x*, the receiver clock offset *dtr*, the zenith wet delay *zwd* and the ionosphere-free ambiguity *λIFNIF*. By adopting a Kalman filter, the unknown parameters can be exactly estimated. It should be noted that in this paper the ionosphericfree Doppler combination is used to derive the velocity [23].

#### *2.2. INS Model*

The mechanization of INS in the navigation frame (*n*-frame) can be expressed as an integral process for the following equation [33,34]:

$$
\begin{bmatrix}
\dot{\boldsymbol{p}}\_{INS}^{\boldsymbol{n}} \\
\dot{\boldsymbol{\omega}}\_{INS}^{\boldsymbol{n}} \\
\dot{\boldsymbol{\mathcal{C}}}\_{b}^{\boldsymbol{n}}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{D}^{-1} \boldsymbol{v}\_{INS}^{\boldsymbol{n}} \\
\boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}} \boldsymbol{f}^{b} - \left(2\boldsymbol{\omega}\_{\dot{ic}}^{\boldsymbol{n}} + \boldsymbol{\omega}\_{c\boldsymbol{n}}^{\boldsymbol{n}}\right) \times \boldsymbol{v}\_{INS}^{\boldsymbol{n}} + \mathbf{g}^{\boldsymbol{n}} \\
\boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}} \left(\boldsymbol{\omega}\_{\dot{ib}}^{\boldsymbol{b}} \times\right) - \left(\boldsymbol{\omega}\_{\dot{ic}}^{\boldsymbol{n}} + \boldsymbol{\omega}\_{c\boldsymbol{n}}^{\boldsymbol{n}}\right) \times \boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}}
\end{bmatrix} \tag{8}
$$

where *D*−1 is the diagonal matrix to transform the rectangular coordinates into geodetic coordinates; *pn INS* and *vn INS* denote the INS position vector and INS velocity vector in *n*-frame; *Cn b* represents the transformation matrix from body-frame (*b*-frame) to *n*-frame; *f b* is the inertial measurement unit (IMU) specific force measurement projected in *b*-frame; *ωb ib* denotes the IMU angular rate measurement expressed in *b*-frame; *ωn ie* is the vector of earth rotation rate expressed in *n*-frame; *gn* is the gravity vector presented in *n*-frame; *ωn en* represents the rotation rate of ECEF frame relative to that of *n*-frame expressed in *n*-frame; The symbol "×" represents the cross-product operator. Obviously, the INS update of position, velocity, and attitude can be implement based on Equation (8).

The perturbation of position, velocity and attitude must be considered in the INS data processing. In this paper, the error model of attitude, velocity, and position is defined as the psi-angle error model, which reads as

$$
\begin{bmatrix}
\delta \dot{p}\_{INS}^{\text{u}} \\
\delta \dot{\boldsymbol{w}}\_{INS}^{\text{u}} \\
\dot{\Phi}
\end{bmatrix} = \begin{bmatrix}
\end{bmatrix} \tag{9}
$$

where *δpn INS* and *δvn INS* are the error vector of INS position and velocity, respectively; *φ* indicates the vector of misalignment angles; *δgn* is the gravity error; *δf b* is error vector of specific force measurement; *δωb ib* denotes error vector of the IMU angular rate. Generally, the error vectors of specific force and angular rate measurement include scale factor errors, bias errors, and white noise [35]. Here, only bias errors are considered and the bias errors of IMU sensors are described as the first Gauss-Markov procedure [34]:

$$
\begin{bmatrix}
\delta \dot{\mathcal{B}}\_a \\
\delta \dot{\mathcal{B}}\_\mathcal{S}
\end{bmatrix} = -\frac{1}{\pi} \begin{bmatrix}
\delta \mathcal{B}\_a \\
\delta \mathcal{B}\_\mathcal{S}
\end{bmatrix} + \begin{bmatrix}
w\_a \\ w\_\mathcal{S}
\end{bmatrix} \tag{10}
$$

where *δBg* and *δBa* represent the error vector of gyro biases and accelerometer biases, respectively; *τ* is the correction time; and *wa*, *wg* represents the driving noise for accelerometer and gyro, respectively.

#### *2.3. PPP-B2b/INS Loosely Coupled Integration Model*

A typical state vector of 15 states is used in PPP-B2b/INS loosely coupled integration, which can be expressed as

$$\mathbf{X}\_k = \begin{bmatrix} \delta p\_{INS}^n & \delta \boldsymbol{\sigma}\_{INS}^n & \boldsymbol{\Phi} & \delta \mathbf{B}\_a & \delta \mathbf{B}\_\circ \end{bmatrix}^T \tag{11}$$

The system model of PPP-B2b/INS loosely coupled integration in the discrete form can be simplified as

$$\mathbf{X}\_{k} = (\mathbf{I} + \mathbf{F} \cdot \mathbf{\hat{\boldsymbol{\omega}}} t) \mathbf{X}\_{k-1} + \mathbf{\hat{\Gamma}}\_{k-1} \mathbf{w}\_{k-1}, \mathbf{w}\_{k-1} \sim N(\mathbf{0}, \mathbf{Q}\_{k-1}) \tag{12}$$

where *Xk*−1 and *Xk* are the vector of state parameter at epoch *k* − 1 and *k*, respectively; *I* is the identify matrix; Δ*t* denotes the time interval between two adjacent epochs; *wk*−1 represents the vector of the process noise at epoch *k* − 1, and **Γ***k*−1 denotes the matrix of the noise distribution; *Qk*−1 represents the covariance matrix of process noise; *F* is the dynamic matrix, which can be derived from the INS error model as follows:

$$F = \begin{bmatrix} F\_{rr} & F\_{rv} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ F\_{vv} & F\_{vv} & \left(\mathbf{f}^b \times \right) & \mathbf{C}\_b^n & \mathbf{0} \\ F\_{\phi r} & F\_{\phi v} & -\left(\omega\_{ic}^n + \omega\_{cn}^n \times \right) & \mathbf{0} & -\mathbf{C}\_b^n \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & -\frac{1}{\tau}I & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & -\frac{1}{\tau}I \end{bmatrix} \tag{13}$$

where **0** denotes the zero matrix; *Frr*, *Frv*, *Fvr*, *Fvv*, *Fφr*, and *Fφv* are sub-matrices of dynamic matrix, which can be derived from Equation (9). Then the Kalman time update can be expressed as [36]

$$\begin{cases} X\_k^- = \Phi\_{k,k-1} X\_{k-1}^+ \\ P\_k^- = \Phi\_{k,k-1} P\_{k-1}^+ \Phi\_{k,k-1}^T + Q\_{k-1} \end{cases} \tag{14}$$

where the superscript "+" and "−" denote the estimated and predicted information; **Φ***k*,*k*−1 = (*I* + *F*·Δ*t*) represents the state transition matrix; *Pk* represents the covariance matrix of the states.

When the PPP-B2b positioning solution and INS predicted solution are available at the same epoch, the Kalman filter measurement update can be operated. However, the positions and velocities from INS are referred to the IMU center, while those from PPP-B2b positioning are typically based on the GNSS antenna phase center. Therefore, the corresponding lever-arm offsets must be compensated. Here, the accurate lever-arm offsets are measured and directly applied to the position and velocity from INS. Then, the observation model for PPP-B2b/INS loosely coupled integration can be expressed as

$$\mathbf{Z}\_k = H\_k \mathbf{X}\_k + \begin{bmatrix} \mathfrak{e}\_p \\ \mathfrak{e}\_v \end{bmatrix}\_{\prime} \begin{bmatrix} \mathfrak{e}\_p \\ \mathfrak{e}\_v \end{bmatrix} \sim N(0, \mathbf{R}\_k) \tag{15}$$

with

$$\mathbf{Z}\_{k} = \begin{bmatrix} \mathbf{p}\_{GNSS}^{\mathrm{n}}\\ \mathbf{v}\_{GNSS}^{\mathrm{n}} \end{bmatrix} - \begin{bmatrix} \mathbf{p}\_{INS}^{\mathrm{n}}\\ \mathbf{v}\_{INS}^{\mathrm{n}} \end{bmatrix} \tag{16}$$

$$H\_k = \begin{bmatrix} I\_{3 \times 3} & \mathbf{0}\_{3 \times 3} & \left(\mathbf{C}\_b^n I^b \times \right) & \mathbf{0}\_{3 \times 3} & \mathbf{0}\_{3 \times 3} \\ \mathbf{0}\_{3 \times 3} & I\_{3 \times 3} & H\_\phi & \mathbf{0}\_{3 \times 3} & -\mathbf{C}\_b^n \left(I^b \times \right) \end{bmatrix} \tag{17}$$

where *Hφ* = *ωn ie* + *ωn en* × *Cn b* - *l b*× + *Cn b* - *l b* × *ωb ib* × ; *pn GNSS*, *vn GNSS* denote the vector of position and velocity obtained from PPP-B2b positioning; *εp εv T* denotes the vector of measurement noise of position and velocity and *Rk* represents its related covariance matrix; *l b* represents the vector of lever-arm offsets; *Hk* is the design matrix; *Zk* is defined as the innovation vector, which consists of the difference of positions and velocities between PPP-B2b positioning solution and INS predicted solution. Finally, the Kalman measurement update can be given as

$$\begin{cases} \mathbf{K}\_k = \mathbf{P}\_k^- \mathbf{H}\_k^T \left( \mathbf{H}\_k \mathbf{P}\_k^- \mathbf{H}\_k^T + \mathbf{R}\_k \right)^{-1} \\ \mathbf{X}\_k^+ = \mathbf{X}\_k^- + \mathbf{K}\_k (\mathbf{Z}\_k - \mathbf{H}\_k \mathbf{X}\_k^-) \\ \mathbf{P}\_k^+ = (\mathbf{I} - \mathbf{K}\_k \mathbf{H}\_k) \mathbf{P}\_k^- \end{cases} \tag{18}$$

where *Kk* is the Kalman gain.

### **3. Description of Experiments**

In order to assess the performance of BDS-3 PPP-B2b/INS loosely coupled integration, two land vehicle experiments in open-sky scenery and urban canyon scenery (hereinafter referred to as Experiment A and Experiment B, respectively) were carried out in Qingdao, China. In these two experiments, a tactical-grade IMU (ISA100C) [37] and a MEMS-grade IMU (ADIS-16505) [38] were equipped to evaluate the impact of IMU grade on BDS-3 PPP-B2b/INS loosely coupled integration. Two GNSS antennas were mounted on the roof of the land vehicle, and connected to a NovAtel PwrPak7 receiver [39] and a FRII-PLUS PPP-B2b receiver (http://www.femtomes.com, accessed on 20 October 2021). The FRII-PLUS receiver was only used to collected PPP-B2b messages, while all GNSS observations used in the evaluation were collected by the PwrPak7 receiver. The sampling rate of GNSS observations was 1 Hz. Figure 1 shows the installation of experimental equipment. Relevant specific information for the IMU sensors is listed in Table 1.

**Figure 1.** Experimental equipment mounted on the vehicle.

**Table 1.** Parameters of the IMU sensors used in the experiments.

During these two experiments, a GNSS base station was established on the roof of Engineering-C building in the China University of Petroleum (East China). It should be noted that the distance between the base and rover stations was less than 5 km. Using GNSS observations collected by PwrPak7, and IMU data from ISA100C, together with GNSS measurements at the base station, the smoothed RTK/INS tightly coupled solution was obtained by commercial Inertial Explorer software (IE 8.90) from the NovAtel company, and this solution was used as reference.

In PPP-B2b positioning, the B1C and B2a signals of BDS-3, L1 and L2 signals of GPS were selected to form the ionospheric-free combination, and the cut-off elevation angle was set to 7 degrees. The remaining processing strategies are shown in detail in Table 2.

**Table 2.** Strategies of PPP-B2b positioning.

#### **4. Result and Discussion**

In this section, the performance of BDS-3 PPP-B2b/INS loosely coupled integration was investigated in both open-sky and complex urban environments.

#### *4.1. Experiment A*

Experiment A was conducted from 16:30:40 to 16:59:59 on 3 December 2021 in GPS time. The trajectory of the land vehicle is shown in Figure 2.

During Experiment A, the velocities of the land vehicle were within ±20 m/s in the east and north directions and within ±1 m/s in the vertical direction. The number of visible BDS-3/GPS satellites and the position dilution of precision (PDOP) values are shown in Figure 3. The number of visible satellites ranged from 9 to 15 with an average value of 12.5, and the PDOP varied from 1.28 to 2.70 with an average value of 1.59.

Figure 4 shows the positioning errors of PPP-B2b only, PPP-B2b/MEMS-IMU loosely coupled integration, and PPP-B2b/tactical-IMU loosely coupled integration, respectively. As with traditional PPP processing, a convergence period is indispensable in PPP-B2b positioning and PPP-B2b/INS loosely coupled integration. In this study, the convergence time was defined for horizontal/vertical positioning accuracy better than 30 cm/60 cm, with such a positioning accuracy for at least 60 continuous epochs. The statistics of the positioning results are summarized in Table 3. It should be noted that positioning results before full convergence were excluded when calculating positioning accuracy. The biases

in the east, north and up directions were 6.7/23.7/26.5 cm for PPP-B2b only, and the corresponding values for PPP-B2b/tactical-IMU and PPP-B2b/MEMS-IMU loosely coupled integration schemes were 6.6/23.3/26.0 cm and 6.2/23.5/25.8 cm, respectively. In terms of the root mean square (RMS) value of positioning errors, the improvement of PPP-B2b/INS loosely coupled integration was not significant compared to PPP-B2b only. This is because Experiment A was carried out within an open-sky environment and thus PPP-B2b on its own was already able to obtain high-precision positioning accuracy.

**Figure 2.** The trajectory of land vehicle in Experiment A.

**Figure 3.** Number of visible BDS-3/GPS satellites and PDOPs in Experiment A.

**Figure 4.** Positioning errors of PPP-B2b only (**top** panel), PPP-B2b/MEMS-IMU loosely coupled integration (**middle** panel), and PPP-B2b/tactical−IMU loosely coupled integration (**bottom** panel), for Experiment A.

**Table 3.** Statistics of the positioning errors (unit: cm).

It is recognized that PPP-B2b only is unable to obtain a high-precision solution during GNSS signal outages, while PPP-B2b/INS loosely coupled integration continues to work by INS mechanization. Four periods of GNSS signal outage were simulated to present the superiority of PPP-B2b/INS loosely coupled integration. The periods of GNSS outage were simulated from the 600th/900th/1200th/1500th epoch, and each period lasted 30 s. The corresponding positioning errors of PPP-B2b/MEMS-IMU and PPP-B2b/tactical-IMU loosely coupled integration are plotted in Figures 5 and 6, respectively. Clearly, the positioning accuracy become progressively worse with increased time of GNSS outages in both PPP-B2b/MEMS-IMU and PPP-B2b/tactical-IMU loosely coupled integration schemes. When the GNSS signal outage time increased to 30 s, the positioning errors in the east, north, and up directions decreased to 300.0 cm, 498.0 cm, and 41.0 cm, respectively, for PPP-B2b/MEMS-IMU loosely coupled integration. The results for the PPP-B2b/tactical-IMU loosely coupled integration scheme dropped to 18.6 cm, 21.8 cm, and 6.1 cm in the east, north, and up components, respectively. Table 4 lists the statistics results, showing that the positioning performance of PPP-B2b/tactical-IMU loosely coupled integration is significantly better than that of PPP-B2b/MEMS-IMU loosely coupled integration.

**Figure 5.** Positioning errors of the PPP-B2b/MEMS−IMU loosely coupled integration during the simulated GNSS signal outages.

**Figure 6.** Positioning errors of the PPP-B2b/tactical−IMU loosely coupled integration during the simulated GNSS signal outages.

**Table 4.** Statistics of positioning errors for PPP-B2b/INS loosely coupled integration during GNSS outages (unit: cm).


**4. EVs as drug delivery vehicles**

Today's medicine is increasingly focused on personalized treatment methods, on targeted therapies that act at the molecular level. One of the concepts aimed at these aspects is that of theranostics, which aims at diagnosis, treatment, and concomitant follow-up of the response by using very specific drug delivery systems [76]. In this sense, EVs are an extremely useful tool for passive diagnosis (especially in neoplastic pathologies, through the ability to identify the tumor type based on the miRNA, mRNA, and mitochondrial RNA profile of EVs) and active (by associating EVs with advanced imaging methods). Thus, numerous platforms based on EVs technology have been developed for theranostic purposes, namely, transition metal-labeled exosomes, nanoparticle-loaded exosomes, bioluminescently labeled exosomes, nanocluster loaded exosomes, metabolically labeled exosomes. The main applications at present are those in the oncology field, but the lack of uniformity of clear classifications, increased immunogenicity of EVs, and the limited number of drugs that can be loaded at EVs highlight the need for future studies for widespread application of these diagnostic and therapeutic tools [77].

EVs are used as transporters for a variety of substances, ranging from small molecules, small interfering RNA, mRNA, and microRNAs to drugs with suboptimal pharmaceutical effects, carrying active constituents through biological barriers [4, 11].

Exosomal transporters present advantages, as they can travel efficiently between cells, smoothly passing their cargo along the cell membrane, keeping it biologically active, and crossing hard-to-penetrate barriers, such as the blood-brain barrier. Important issues regarding exosome-based drug delivery vehicles are the precise method of exosome loading, without altering the biological characteristics, and the scalable repeatable production of exosome categories [18, 78].

Exosomes tend to have special homing targets, influenced by their cell of origin [18]. Their membrane can be modified, to amplify the targeting of specific cells [18, 79].

The content of EVs can be loaded exogenously (integration of small proteins, RNA, or other molecules) or endogenously (assuring that cells possess the ways to integrate small molecules, proteins or RNA into EVs during their formation) [4, 80]. Exogenous changes of EVs can be done after their collection, incorporating the cargo into EVs through different methods: coincubation (with no modification of vesicle size distribution or integrity, electroporation, and sonication) [79, 80]. The endogenous loading can be obtained through artificial adjustment of the parental cell to overexpress certain proteins or RNA, that can be integrated into secreted EVs afterward [81].

Human MSCs, multipotent adult progenitors, could be an adequate source of exosomes for drug delivery. Their transplantation has been investigated in numerous trials and proved to be safe, MSCs also produce immunologically inert exosomes [18]. As a delivery system, studies have shown that EV derived from MSCs can transfer therapeutic drugs to diseased cells [19]. Their efficiency is based on the adhesion proteins on their surface (like integrins, extracellular matrix proteins, tetraspanins), which facilitate the penetration of the cellular membrane and the accumulation of EVs in the diseased cells [82]. Other characteristics that make EVs an ideal candidate for a drug delivery system are their decreased toxicity and immunogenicity, as well as their potential to cross the blood-brain barrier [83, 84]. A study conducted by S. Kamerkar *et al*. has shown that EVs produced by MSCs can transfer small interfering RNA targeting the oncogenic KRas (G12D) mutants to pancreatic cancer cells, increasing cells apoptosis and decreasing the risk of metastatic disease [85].

Exosomes have significant immune properties, modulating immunological responses and facilitating antigen presentation [11, 86]. Exosomes derived from dendritic cells can conduct MHC class I/peptide complexes to other dendritic cells for *in vivo* activation of cytotoxic T lymphocytes and promote T cell-dependent antitumor responses *in vivo* [86, 87]. Dendritic cell exosomes have been previously loaded with antigenic peptides, to activate T cell proliferation, with possible use as vaccines against infectious or neoplastic diseases. Due to the immunogenic nature of dendritic cell exosomes, their use as drug delivery vehicles is not ideal. A more suitable choice would be human ESCs-derived mesenchymal cells [86, 88, 89].

Clinical trials with therapies based on EVs are studied in malignancies, such as melanoma, non-small cell lung cancer, colon cancer, metastatic pancreatic cancer, bronchopulmonary dysplasia, malignant ascites, and pleural effusion, but also chronic kidney disease, type 1 diabetes, insulin resistance and chronic inflammation polycystic ovary syndrome, ulcers, and acute ischemic stroke [4]. Exosomes were shown to transport curcumin and chemotherapeutics, such as doxorubicin and paclitaxel [19, 90].

#### *Extracellular Vesicles as Intercellular Communication Vehicles in Regenerative Medicine DOI: http://dx.doi.org/10.5772/intechopen.101530*

The implication of MSCs was also evaluated in patients with anthracyclineinduced cardiomyopathy. Mitochondrial transfer, mediated by large EVs, diminished injury determined by doxorubicin, in patient-specific induced pluripotent SC-derived cardiomyocytes. MSCs could ameliorate cardiac function in anthracycline-induced cardiomyopathy, regardless of regeneration effects [91].

Liposomes possess many favorable characteristics as drug delivery vehicles, being used in the transportation of anti-cancer drugs, anti-fungal medication, and analgesics [92–98]. Liposomes have a phospholipid membrane that helps with the incorporation of hydrophilic or hydrophobic drugs, and they can also deliver the carried drugs to the targeted points through plasma membrane breaching. To diminish the recognition by opsonins and their clearance, liposomes can be covered with polymers (PEG). Their membranes can be adapted, to present ligands or antibody elements, which can interact with specific cells and amplify the targeted drug delivery. Liposomes, with easy-to-control properties, can be loaded with drugs, DNA, diagnostic instruments, enzymes, or peptides. Drugs included in liposomes have attenuated toxicity and do not provoke unwanted toxic reactions [99]. Liposomal drugs have various routes of administration, such as parenteral, oral, topical, and even through aerosols [99].

Synthetic liposomes, although very useful, are overcome by EVs (naturally derived liposomes), which have lower toxicity [19]. Exosomes are considered superior drug delivery vehicles, as an alternative to liposomes. In contrast to the latter, exosomes are usually adequately tolerated by the human body and do not present intrinsic toxicity. They can deliver their content through the plasmatic membrane and protect against its early transformation and elimination [19]. Since exosomes can be found in a variety of biological fluids, such as blood, urine, breast milk the delivered drugs will be well tolerated, less toxic, and with a longer circulating half-life [19, 100].


**4.3 Indonesia financial link with Argentina and Turkey**

## *4.3.1 Measuring contagion with correlation*

The figure above shows that Jakarta composite index (JCI) has the highest correlation with MERVAL, namely, the Argentina exchange (0.9), and the second highest is with ISTANBUL (0.89) or Turkish stock exchange. Thus, the government of Indonesia has to pay attention to the high correlation potential with both countries. Correlation with NASDAQ (US exchange) is significant but slightly below Turkey (0.88). NASDAQ has a high and significant correlation with the three emerging market stock exchanges (Turkey, Indonesia, and Argentina) (**Figure 5** and **Table 4**).

**Table 5** shows the adjusted correlation coefficients before and during the crisis. All of the coefficients increase during the crisis, except the contagion between Indonesia and Malaysia which the stock price fell sharply. Meanwhile, the correlation between the Thai and Indonesian exchange rates and between the Thai and Malaysian exchange rates experiences a high surge during the crisis.

#### **Table 4.**

*Correlation result.*

statistically significant. In fact, there is a statistically significant decrease in the correlation coefficient between the Malaysian and Indonesian stock markets. One plausible explanation for this decline is the adoption of Malaysian capital controls in

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial…*

*BORSA (Turkey)—JCI (Indonesia) correlation 2005–2019. Source: Authors, 2019.*

in some years. Changes in the correlation between the stock market and the

*MERVAL (Argentina)—JCI (Indonesia) correlation 2005–2019. Source: Authors, 2019.*

**Table 6** shows the correlation coefficients of the stock market and exchange rate

exchange rate of Indonesia from 2017 to 2018 were not as big as when the economic crisis of 1997. This happens because after the crises, Indonesia becomes more vigilant so that various agencies were formed to predict and overcome crises. There is an increase in the correlation between Indonesian stock prices from 2017 to 2018 both with Argentina and Turkey's stock prices. Something similar also happens to the exchange rate. There is an increase in the correlation between the Indonesian exchange rate from 2017 to 2018 both with the exchange rates of Argentina and

late August and early September 1998 [24].

*DOI: http://dx.doi.org/10.5772/intechopen.92275*

Turkey (**Table 6**).

**Figure 7.**

**51**

**Figure 6.**

#### **Table 5.**

*Contagion in the currency and stock markets (adjusted correlation coefficients) in 1997 from Iriana and Sjöholm [24].*

A negative and statistically significant coefficient means that there is evidence of contagion. The increases in the correlation between the Thai and Indonesian exchange rate and between the Thai and Malaysian exchange rate are statistically significant. Therefore, the results show that difficulties in Thailand are transmitted to the Indonesian and Malaysian currency markets. There are no signs of contagion on the stock market because an increase in the correlation coefficient is not

#### **Table 6.**

*Potential contagion: correlation rises during crisis.*

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial… DOI: http://dx.doi.org/10.5772/intechopen.92275*

**Figure 6.** *BORSA (Turkey)—JCI (Indonesia) correlation 2005–2019. Source: Authors, 2019.*

statistically significant. In fact, there is a statistically significant decrease in the correlation coefficient between the Malaysian and Indonesian stock markets. One plausible explanation for this decline is the adoption of Malaysian capital controls in late August and early September 1998 [24].

**Table 6** shows the correlation coefficients of the stock market and exchange rate in some years. Changes in the correlation between the stock market and the exchange rate of Indonesia from 2017 to 2018 were not as big as when the economic crisis of 1997. This happens because after the crises, Indonesia becomes more vigilant so that various agencies were formed to predict and overcome crises. There is an increase in the correlation between Indonesian stock prices from 2017 to 2018 both with Argentina and Turkey's stock prices. Something similar also happens to the exchange rate. There is an increase in the correlation between the Indonesian exchange rate from 2017 to 2018 both with the exchange rates of Argentina and Turkey (**Table 6**).

**Figure 7.** *MERVAL (Argentina)—JCI (Indonesia) correlation 2005–2019. Source: Authors, 2019.*

A negative and statistically significant coefficient means that there is evidence of

**Stock market Exchange rate Year JCI-MERVAL JCI-Turkey Year IDR-ARS IDR-TRY** 0.338 0.779 2015 0.341 0.930 0.924 0.154 2016 0.252 0.279 0.917 0.932 2017 0.713 0.725 0.952 0.950 2018 0.841 0.826

**Correlation NASDAQ MERVAL JCI ISTANBUL BROAD\_USD**

12.578 13.295 12.280 — BROAD\_USD 0.025 0.062 0.094 0.208 1.000

**Period Stock market Exchange rate**

**Thailand-Malaysia**

Tranquil 0.25 0.26 0.08 0.02 0.12 0.14 Crisis 0.29 0.07 0.09 0.19 0.14 0.31

**Malaysia-Indonesia**

0.155 0.387 0.593 1.330 —

0.54 2.55\*\* 0.14 2.26\*\* 0.037 2.29\*\*

**Thailand-Indonesia**

**Malaysia-Indonesia** **Thailand-Malaysia**

NASDAQ 1.000

*Public Sector Crisis Management*

*Source: Authors Calculation, 2019.*

**Thailand-Indonesia**

**Table 4.** *Correlation result.*

> Statistically significant

**Table 5.**

*Sjöholm [24].*

*\*\*Significant at a 5% level. Source: Iriana and Sjöholm [24].*

*Source: Authors Calculation, 2019.*

*Potential contagion: correlation rises during crisis.*

**Table 6.**

**50**

— MERVAL 0.964 1.000

22.783 — JCI 0.878 0.901 1.000

11.436 12.995 — ISTANBUL 0.896 0.905 0.891 1.000

contagion. The increases in the correlation between the Thai and Indonesian exchange rate and between the Thai and Malaysian exchange rate are statistically significant. Therefore, the results show that difficulties in Thailand are transmitted to the Indonesian and Malaysian currency markets. There are no signs of contagion

*Contagion in the currency and stock markets (adjusted correlation coefficients) in 1997 from Iriana and*

on the stock market because an increase in the correlation coefficient is not

#### *<IP_ADDRESS> Dynamic conditional correlation*

### *<IP_ADDRESS>.1 Stock market*

Throughout 2005–2019, the correlation between the Turkish stock market and Indonesian stock market is in the range of 0.6–0.9. The biggest decline occurs in 2007 and 2010 which reach 0.3, and the highest correlation increase occurs in 2018 (**Figure 6**). Meanwhile, **Figure 7** shows a correlation between Argentine stock market and Indonesian stock market which experience an upward trend since 2011. The upward trend between Argentina's stock market and Indonesia's stock market also happened since 2013, as well as in the periods of 2016–2018.

#### *<IP_ADDRESS>.2 Exchange rate*

The exchange rate correlations between 2005 and 2019 experience fluctuations. The correlation between lira and rupiah begins with a negative value, which is almost touching 0.5. In addition to 2005, in 2008 and 2014 lira and rupiah are also

negatively correlated. Meanwhile, the correlation between peso and rupiah experiences a negative correlation from the end of 2008 to the beginning of 2009 and at

**Correlation Argentina Indonesia Turkey 2002Q3–2008Q4** Argentina 1.000 0.036 0.125 Indonesia **0.036** 1.000 0.046 Turkey 0.125 0.046 1.000 **2009Q1 - 2018Q4** Argentina 1.000 0.207 0.088 Indonesia 0.207 1.000 0.423 Turkey 0.088 **0.423** 1.000

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial…*

**Table 7** shows the net correlation after controlling fundamental between Argentina, Indonesia, and Turkey in terms of capital markets. The results show that the highest correlation was between Indonesia and Turkey from 2009 to 2018. Meanwhile, the lowest correlation was between Indonesia and Argentina in 2002

By using monthly data and the lira exchange rate as the endogenous variables, the research will analyze the impact of weaker lira on the Indonesian economy (VAR dynamic method). **Table 8** shows that the change in lira exchange rate is

**Rupiah Interbank**

**rates**

**Net export** **LIP (production index)**

dominant in explaining changes in the rupiah exchange rate (26–37%).

1 0.031 26.635 0.096 73.269 0.000 0.000 0.000 2 0.033 35.264 0.080 63.237 0.177 0.852 0.390 3 0.033 35.947 0.216 62.291 0.174 0.895 0.476 4 0.034 37.063 0.224 61.157 0.174 0.881 0.500 5 0.034 37.077 0.263 61.073 0.174 0.909 0.504 6 0.034 37.088 0.264 61.058 0.175 0.910 0.506 7 0.034 37.092 0.266 61.049 0.175 0.913 0.505 8 0.034 37.092 0.266 61.048 0.175 0.914 0.506 9 0.034 37.092 0.266 61.048 0.175 0.914 0.506 10 0.034 37.092 0.266 61.048 0.175 0.914 0.506

**purchase**

the beginning of 2013 (**Figures 8** and **9**).

*DOI: http://dx.doi.org/10.5772/intechopen.92275*

*Net correlation in stock market after extracting fundamental.*

quarter 3 to the end of 2008.

**Period SE Lira Net foreign**

*Variance decomposition of rupiah exchange rate (level).*

*Source: Author Calculation, 2019.*

**Table 7.**

*4.3.2 VAR results*

*Source: Authors, 2019.*

**Table 8.**

**53**

**Figure 8.** *Turkey lira—Indonesia rupiah correlation 2005–2019. Source: Authors, 2019.*

**Figure 9.** *Peso Argentina—Indonesia rupiah correlation 2005–2019. Source: Authors, 2019.*

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial… DOI: http://dx.doi.org/10.5772/intechopen.92275*

### **Table 7.**

*<IP_ADDRESS> Dynamic conditional correlation*

*Public Sector Crisis Management*

Throughout 2005–2019, the correlation between the Turkish stock market and Indonesian stock market is in the range of 0.6–0.9. The biggest decline occurs in 2007 and 2010 which reach 0.3, and the highest correlation increase occurs in 2018 (**Figure 6**). Meanwhile, **Figure 7** shows a correlation between Argentine stock market and Indonesian stock market which experience an upward trend since 2011.

The exchange rate correlations between 2005 and 2019 experience fluctuations.

The correlation between lira and rupiah begins with a negative value, which is almost touching 0.5. In addition to 2005, in 2008 and 2014 lira and rupiah are also

*Turkey lira—Indonesia rupiah correlation 2005–2019. Source: Authors, 2019.*

*Peso Argentina—Indonesia rupiah correlation 2005–2019. Source: Authors, 2019.*

The upward trend between Argentina's stock market and Indonesia's stock market also happened since 2013, as well as in the periods of 2016–2018.

*<IP_ADDRESS>.1 Stock market*

*<IP_ADDRESS>.2 Exchange rate*

**Figure 8.**

**Figure 9.**

**52**

*Net correlation in stock market after extracting fundamental.*

negatively correlated. Meanwhile, the correlation between peso and rupiah experiences a negative correlation from the end of 2008 to the beginning of 2009 and at the beginning of 2013 (**Figures 8** and **9**).

**Table 7** shows the net correlation after controlling fundamental between Argentina, Indonesia, and Turkey in terms of capital markets. The results show that the highest correlation was between Indonesia and Turkey from 2009 to 2018. Meanwhile, the lowest correlation was between Indonesia and Argentina in 2002 quarter 3 to the end of 2008.

#### *4.3.2 VAR results*

By using monthly data and the lira exchange rate as the endogenous variables, the research will analyze the impact of weaker lira on the Indonesian economy (VAR dynamic method). **Table 8** shows that the change in lira exchange rate is dominant in explaining changes in the rupiah exchange rate (26–37%).

#### **Table 8.**

*Variance decomposition of rupiah exchange rate (level).*

### *Public Sector Crisis Management*

and NFP. The decline of the lira may increase in the interbank rates, but there is a decline in NFP at the start, even though the NFP could recover after the third month. The weakening of lira also caused a decline in net exports and a manufacturing production index even though with little impact (**Figure 10**).

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial…*

sensitive to get the impact of capital flow.

*DOI: http://dx.doi.org/10.5772/intechopen.92275*

*R-square decomposition based on Shepley decomposition. Source: [27].*

**4.4 Rupiah exchange rate model results**

**4.5 Exchange rate model result**

dummy caused rupiah to depreciate) (**Table 11**).

þ *ϵt*

(**Figure 12**).

**55**

**Figure 11.**

Claessens [27] explains that the main reasons of fragility are liquidity and investor based, while macrofundamentals only have a little explaining power except for bond (**Figure 11**). According to **Table 10**, Indonesia, Turkey, and Argentina are

Indonesian capital and financial account in the fourth quarter of 2018 showed a good performance and even had obtained the highest value since 2012. Unfortunately, this is not supported by the current account condition. The current account in the fourth quarter of 2018 continued to deteriorate compared to the previous periods. This occurs because of the global economic slowdown that is currently happening so that Indonesia's exports to several countries have decreased

*EXCH*\_*RATEt* ¼ *α*10 þ *β*11*TBt*�1 þ *β*12*FINANCIALt*�1 þ *β*13*PRIMARY INCOMEt*�1 þ *β*14*LIRAt* þ *β*15*DUMMYt* þ *β*16*AR*ð Þ1 *t* þ *ϵt*

Based on the regression results, it can be concluded that Indonesia's trade balance and financial account have a significant impact on the exchange rate on α 0.1 with a negative coefficient (the surplus of financial account caused rupiah to appreciate). In contrast, lira and crisis dummy have a significant positive effect on exchange rates (Lira depreciation followed by rupiah depreciation, and crisis

*LOG EXCH* ð Þ \_*RATE t* ¼ *α*10 þ *β*11*TBt* þ *β*12*LOG M*ð Þ2 *t* þ *β*13*DUMMYt* þ *β*14*AR*ð Þ1 *t*

Based on **Table 12**, it can be concluded that broad money (M2) and dummy variables significantly influence the exchange rate. Thus, if there is a higher amount

(16)

(17)

**Figure 10.** *IRF graph. Source: Authors, 2019.*

#### **Table 9.**

*Variance decomposition of Indonesia's industrial production index.*

The changes in the Indonesia production index in response to changes in the lira exchange rate variable are relatively small. Interbank rates, net exports, and rupiah exchange rates (apart from the production index itself) get the highest change. So we can conclude that the contagion for Turkey has hit Indonesia more on the financial market, especially exchange rate, and has a small effect on real sector activity represented by industrial production index (**Figure 10** and **Table 9**).

The weakening lira is proven to be the reason why the rupiah moves downward, but the shock would disappear after 5 months (the rupiah value returned to balance). Not only rupiah but also weaker lira also has an impact on the interbank rates *Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial… DOI: http://dx.doi.org/10.5772/intechopen.92275*

**Figure 11.**

*R-square decomposition based on Shepley decomposition. Source: [27].*

and NFP. The decline of the lira may increase in the interbank rates, but there is a decline in NFP at the start, even though the NFP could recover after the third month. The weakening of lira also caused a decline in net exports and a manufacturing production index even though with little impact (**Figure 10**).

Claessens [27] explains that the main reasons of fragility are liquidity and investor based, while macrofundamentals only have a little explaining power except for bond (**Figure 11**). According to **Table 10**, Indonesia, Turkey, and Argentina are sensitive to get the impact of capital flow.

#### **4.4 Rupiah exchange rate model results**

Indonesian capital and financial account in the fourth quarter of 2018 showed a good performance and even had obtained the highest value since 2012. Unfortunately, this is not supported by the current account condition. The current account in the fourth quarter of 2018 continued to deteriorate compared to the previous periods. This occurs because of the global economic slowdown that is currently happening so that Indonesia's exports to several countries have decreased (**Figure 12**).

#### **4.5 Exchange rate model result**

$$\begin{aligned} \text{EXCH\\_RATE}\_t &= a\_{10} + \beta\_{11} \text{TB}\_{t-1} + \beta\_{12} \text{FINANCE}\_{t-1} + \beta\_{13} \text{PRIMARY INCOMP}\_{t-1} \\ &+ \beta\_{14} \text{LIRA}\_t + \beta\_{15} \text{DUMMWY}\_t + \beta\_{16} \text{AR}(\mathbf{1})\_t + \varepsilon\_t \end{aligned} \tag{16}$$

Based on the regression results, it can be concluded that Indonesia's trade balance and financial account have a significant impact on the exchange rate on α 0.1 with a negative coefficient (the surplus of financial account caused rupiah to appreciate). In contrast, lira and crisis dummy have a significant positive effect on exchange rates (Lira depreciation followed by rupiah depreciation, and crisis dummy caused rupiah to depreciate) (**Table 11**).

$$\begin{aligned} \text{LOG}(\text{EXCH\\_RATE})\_t &= a\_{10} + \beta\_{11} \text{TB}\_l + \beta\_{12} \text{LOG}(\text{M\!2})\_t + \beta\_{13} \text{LUMM} \text{Y}\_t + \beta\_{14} \text{AR}(\text{1})\_t \\ &+ a\_l \end{aligned} \tag{17}$$

Based on **Table 12**, it can be concluded that broad money (M2) and dummy variables significantly influence the exchange rate. Thus, if there is a higher amount

The changes in the Indonesia production index in response to changes in the lira exchange rate variable are relatively small. Interbank rates, net exports, and rupiah exchange rates (apart from the production index itself) get the highest change. So we can conclude that the contagion for Turkey has hit Indonesia more on the financial market, especially exchange rate, and has a small effect on real sector activity represented by industrial production index (**Figure 10** and **Table 9**).

**Figure 10.**

*IRF graph. Source: Authors, 2019.*

*Public Sector Crisis Management*

**Period SE D**

*Source: Authors Calculation, 2019.*

**Table 9.**

**54**

**(LER\_LIRA)**

*Variance decomposition of Indonesia's industrial production index.*

**D (NFP)**

**D (LER)**

1 6080.729 0.222 1.287 0.152 1.396 3.926 93.015 2 6982.743 0.202 1.954 0.349 1.376 4.916 91.202 3 7130.040 0.223 1.866 4.273 3.866 4.797 84.973 4 7163.440 0.672 1.781 4.370 6.371 5.101 81.703 5 7169.733 0.721 1.772 4.351 6.560 5.234 81.363 6 7173.949 0.802 1.775 4.366 6.739 5.223 81.095 7 7176.103 0.811 1.779 4.367 7.014 5.250 80.778 8 7176.375 0.812 1.780 4.366 7.051 5.260 80.731 9 7176.617 0.814 1.780 4.366 7.060 5.260 80.719 10 7176.858 0.815 1.782 4.366 7.085 5.263 80.688

**D (PUAB)**

**D (NET\_EXPORT)**

**D (LIP\_INA)**

The weakening lira is proven to be the reason why the rupiah moves downward, but the shock would disappear after 5 months (the rupiah value returned to balance). Not only rupiah but also weaker lira also has an impact on the interbank rates

#### **Table 10.**

*Results of variance decompositions.*

#### **Figure 12.**

*Current account and financial account 2012Q1–2018Q4. Source: Economics and Finance Statistics, Central Bank of Indonesia.*

of money circulating in the community, then there is a decline in the rupiah exchange rate.

$$\begin{split} \text{LOG}(\text{EXCH\\_RATE})\_t &= a\_{10} + \beta\_{11} \text{FINANCE}L\_t + \beta\_{12} \text{COMPICE}\_t + + \beta\_{13} \text{LOG}(\text{M\!2})\_t \\ &+ \beta\_{14} \text{DUMMY}\_t + \beta\_{15} \text{AR}(\text{1})\_t + \varepsilon\_t \end{split}$$

$$\bf{(18)}$$

Based on the table above, it can be concluded that the commodity price variable

**Table 14** shows that financial accounts significantly influence the exchange rate at a significance level of 10% with a negative direction. The increase in financial accounts surplus can cause rupiah appreciation. The condition of the peso exchange

þ *β*14ARS*t* þ *β*15*DUMMYt* þ *β*16*AR*ð Þ1 *t* þ *ϵt* (19)

is significant at a significance level of 5%. The impact of commodity price to exchange rate is negative; thus when there is an increase in commodity prices by 1%, the value of the rupiah will appreciate by 0.876%. On the other hand, the broad money and crisis dummy have a significant and positive coefficient, which means that if money supply increases and crisis happens, rupiah exchange rate will depre-

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial…*

*EXCH*\_*RATEt* ¼ *α*10 þ *β*11*TBt* þ *β*12*FINANCIALt* þ *β*13*PRIMARYINCOMEt*�2

**Variable Coef. Prob** C 8450.149 0.000 TB (�1) �0.052 0.078\* FINANCIAL (�1) �0.016 0.095\* PRIMARY\_INCOME (�1) �0.071 0.256 LIRA 9.022.642 0.000\*\*\* DUMMY 1,084,258 0.003\*\*\* AR (1) 0.872 0.000\*\*\* Adj R2 0.968 Prob (F-stat) 0.000

*Exchange rate determination regression results – Full model with lira variables.*

*Exchange rate determination regression results – Trade balance only.*

**Variable Coef. Prob** C 7937.000 0.000 TB 0.000 0.103 LOG(M2) 0.089 0.061\* DUMMY 0.155 0.002\*\*\* AR (1) 0.787 0.000\*\*\* Adj R2 0.917 Prob (F-stat) 0.000

ciate (**Table 13**).

*DOI: http://dx.doi.org/10.5772/intechopen.92275*

*\**

*\**

**57**

**Table 12.**

*Significant at α 0.1. \*\*Significant at α 0.05. \*\*\*Significant at α 0.01. Source: Authors, 2019.*

**Table 11.**

*Significant at α 0.1. \*\*Significant at α 0.05. \*\*\*Significant at α 0.01. Source: Authors, 2019.*

*Contagion, Exchange Rate, and Financial Volatility: Indonesian Case in Global Financial… DOI: http://dx.doi.org/10.5772/intechopen.92275*

Based on the table above, it can be concluded that the commodity price variable is significant at a significance level of 5%. The impact of commodity price to exchange rate is negative; thus when there is an increase in commodity prices by 1%, the value of the rupiah will appreciate by 0.876%. On the other hand, the broad money and crisis dummy have a significant and positive coefficient, which means that if money supply increases and crisis happens, rupiah exchange rate will depreciate (**Table 13**).

$$\begin{array}{c} \text{EXCH\\_RATE}\_t = a\_{10} + \beta\_{11} \text{TB}\_t + \beta\_{12} \text{FINANCLAL}\_t + \beta\_{13} \text{PRIAMRY}\_{\text{INCOME}\_t - 2} \\ \quad + \beta\_{14} \text{ARS}\_t + \beta\_{15} \text{DUMMY}\_t + \beta\_{16} \text{AR}(\mathbf{1})\_t + \epsilon\_t \end{array} \tag{19}$$

**Table 14** shows that financial accounts significantly influence the exchange rate at a significance level of 10% with a negative direction. The increase in financial accounts surplus can cause rupiah appreciation. The condition of the peso exchange

#### **Table 11.**

*Exchange rate determination regression results – Full model with lira variables.*

#### **Table 12.**

*Exchange rate determination regression results – Trade balance only.*

of money circulating in the community, then there is a decline in the rupiah

*LOG EXCH* ð Þ \_*RATE t* ¼ *α*10 þ *β*11*FINANCIALt* þ *β*12*COMPRICEt* þ þ*β*13*LOG M*ð Þ2 *t* þ *β*14*DUMMYt* þ *β*15*AR*ð Þ1 *t* þ *ϵt*

*Current account and financial account 2012Q1–2018Q4. Source: Economics and Finance Statistics, Central*

**Country Equity Bond Bank** Turkey 0.56 0.42 0.42 Argentina 0.37 0.14 0.32 Indonesia 0.51 0.69 0.43 South Africa 0.46 0.58 0.50 Israel 0.17 0.36 �0.03 Brazil 0.58 0.52 0.46 Chile �0.06 0.15 0.19 Colombia 0.16 0.02 0.23 Mexico 0.30 0.38 0.27 Peru 0.27 0.33 0.45 Uruguay �0.09 0.44 0.02 Venezuela, Rep. Bol. �0.06 0.29 �0.18 India 0.67 0.16 0.23 China PR: Mainland 0.41 �0.08 0.57 Korea 0.49 0.27 0.43 Malaysia 0.38 0.29 0.45 Pakistan 0.90 0.40 0.12 Philippines 0.64 0.36 0.19 Thailand 0.58 0.36 0.40

(18)

exchange rate.

*Bank of Indonesia.*

**Figure 12.**

**56**

*Source: Claessens [27].*

*Results of variance decompositions.*

*Public Sector Crisis Management*

**Table 10.**

*\*\*\*Significant at α 0.01. Source: Authors, 2019.*


**6. Prosthetic and facial reconstruction**

Oral and facial deformity can cause functional and psychological deterioration in oncological patients. The aims of prosthodontic reconstruction are the rehabilitation of the shape, function and esthetics of the lost anatomical structures by means of artificial substitutes. The main facial subunits that require reconstruction due to malignant pathological processes involve the ear, forehead, eyes and brow, nose, check, lips and chin [33].

Currently, there are multiple workflows where they combine conventional prosthetic preparation with the use of facial scanners and custom 3D impressions

**169**

**Figure 6.**

*A Review of Maxillofacial Rehabilitation Using Osseointegrated Implants in Oncological…*

that allow the direct or indirect manufacture of definitive maxillofacial prostheses

In early clinical stages, surgery is the first decision. However, adjuvant radio-

involved (<1 mm), and in suspicion or confirmation of lymph node metastasis with or without extracapsular extension [36]. In general, radiation therapy includes conventional radiation therapy or intensity modulated radiation therapy (IMRT). The latter is more convenient by precisely targeting radiation to a specific area and reducing the dose to nearby anatomical structures such as malar bone, grafts,

Schiegnitz's study suggests that radiation negatively affects implant survival, but there is no statistically significant difference in survival when implants are inserted 12 months before or after radiation therapy [38]. Other authors have shown favorable results. With implantation at least 6 weeks before radiotherapy since there is a surgical area with less hypoxia, hypovascularity and hypocellularity [27, 39] (**Figure 6**).

Systemic antineoplastic therapy includes induction, neoadjuvant, adjuvant, and palliative therapy [40]. The survival rate of osseointegrated implants in the patient with antineoplastic and antiangiogenic therapy has not yet been elucidated. Controlled animal studies have shown negative effects on the osseointegration process with chemotherapeutic regimens with cisplatin, bevacizumab and sunitinib [41–43]. However, retrospective studies in humans have reported that chemotherapy has no detrimental effects on the osseointegration and functional stability of dental implants.

The success rates reported by these studies were 97.6 and 99.1% [44–46].

therapy is sometimes indicated in cases of close excision margins (<5 mm),

The main advantage with dental prosthodontic reconstruction is the adaptation of a retained implant obturator that allows a posterior palatal seal. With this, most patients successfully recover phonation, swallowing, and chewing function. Artificial facial prostheses allow us to better characterize the demanding anatomy of the face. In most cases, the retention of facial prostheses occurs when copying the anatomical defect, on other occasions magnetized attachments can be used in prosthetic reconstruction or protective glasses that also improve the characterization of the face (**Figure 5**)**.**

*DOI: http://dx.doi.org/10.5772/intechopen.93224*

**7. Radiotherapy treatment and implants**

implants, salivary glands, eyes, and spinal cord [37].

**8. Chemotherapy treatment and implants**

*Histogram showing trajectories, dose and volume used for adjuvant radiotherapy.*

from conventional silicone [34, 35].

**Figure 5.** *Photograph showing implant-supported obturator placed and final photograph.*

*A Review of Maxillofacial Rehabilitation Using Osseointegrated Implants in Oncological… DOI: http://dx.doi.org/10.5772/intechopen.93224*

that allow the direct or indirect manufacture of definitive maxillofacial prostheses from conventional silicone [34, 35].

The main advantage with dental prosthodontic reconstruction is the adaptation of a retained implant obturator that allows a posterior palatal seal. With this, most patients successfully recover phonation, swallowing, and chewing function. Artificial facial prostheses allow us to better characterize the demanding anatomy of the face. In most cases, the retention of facial prostheses occurs when copying the anatomical defect, on other occasions magnetized attachments can be used in prosthetic reconstruction or protective glasses that also improve the characterization of the face (**Figure 5**)**.**


*Article* **Graph, Spectra, Control and Epidemics: An Example with a SEIR Model**

**Giacomo Aletti †, Alessandro Benfenati \*,† and Giovanni Naldi †**

Environmental Science and Policy Department, Università degli Studi di Milano, 20133 Milan, Italy<EMAIL_ADDRESS>(G.A<EMAIL_ADDRESS>(G.N.)

**\*** Correspondence<EMAIL_ADDRESS>
† The three authors are members of the Italian Group GNCS of the Italian Institute "Istituto Nazionale di Alta Matematica", and of the ADAMSS Center of the Università degli Studi di Milano (Italy).

**Abstract:** Networks and graphs offer a suitable and powerful framework for studying the spread of infection in human and animal populations. In the case of a heterogeneous population, the social contact network has a pivotal role in the analysis of directly transmitted infectious diseases. The literature presents several works where network-based models encompass realistic features (such as contacts networks or host–pathogen biological data), but analytical results are nonetheless scarce. As a significant example, in this paper, we develop a multi-group version of the epidemiological SEIR population-based model. Each group can represent a social subpopulation with the same habits or a group of geographically localized people. We consider also heterogeneity in the weighting of contacts between two groups. As a simple application, we propose a simple control algorithm in which we optimize the connection weights in order to minimize the combination between an economic cost and a social cost. Some numerical simulations are also provided.

**Keywords:** epidemic spread; multi-group models; network based model; control of spread dynamics

#### **1. Introduction**

The epidemiological modeling of infectious disease transmission has a long history in mathematical biology, for humans [1–7], animals [3,8] and plants [9–11]. In recent years it has had an increasing influence on the theory and practice of disease management and control, e.g., [12–18]. Indeed the forecast of the spread of an infectious disease is critical to public health decision making.

The proper modeling and analysis of the dynamics of infectious diseases has been a long-standing area of research among many different fields, including economics, social sciences, mathematical biology, physics, computer science and engineering [5]. In the classical population approach, the underlying common factor is the partitioning of the population into "compartments"; we assume that the populations in the various compartments are homogenous in the sense that all individuals behave similarly. The two most common compartments that exist in almost all epidemic models are susceptible (S) and infected (I) [2,3]. The subpopulation S represents individuals who are healthy but susceptible to becoming infected, while I represents individuals who became infected but are able to recover. If the model contains only these two compartments, a given population is initially divided into them. From this basic compartmentalization, there are numerous ways for introducing different interactions within the population. Most of these models for the disease evolution make two basic assumptions. The first assumption states that the population is well-mixed. In such a population, each individual has the same probability of encountering other infected individuals, and thus the resulting force of infection is equal for all. The second assumption states that there are a priori constraints upon the biological process, whilst gradual but random mutation of disease traits (such as transmission rate and infectious period) could occur. More refined epidemic models are required; the entire

**Citation:** Aletti, G.; Benfenati, A.; Naldi, G. Graph, Spectra, Control and Epidemics: An Example with a SEIR Model. *Mathematics* **2021**, *9*, 2987. https://doi.org/10.3390/math9222987

Academic Editors: Mihaela Neamt,u, Eva Kaslik and Anca R ˘adulescu

Received: 31 October 2021 Accepted: 19 November 2021 Published: 22 November 2021

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

**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/).

population cannot just be divided into two or more compartments/groups which are defined by a single quantity. In this paper, we consider the effect of the heterogeneity in the weighting of contacts between two individuals. Moreover, we focus on a meta-population model where the population is previously subdivided into subpopulations that can consist in spatially distinct groups of individuals (neighborhoods, towns, cities, etc.) or groups of individuals with different features. The resulting model is described by a dynamic system defined on a network (graph). We have also added the possibility of varying the weight of the connections between groups in order to formalize the problem of controlling the spread of the epidemic on the network. This generality could also allow the changing of disease features.

The paper is organized as follows: Section 2 introduces the model and the analysis of the corresponding dynamical system. Moreover, as an application of this approach, we will also discuss a new definition of control problem about the spread of epidemics on the network. Section 3 is devoted to numerical experiments using a reduced network with different features. Following the findings of the case study and of previous analysis, the conclusions are presented in Section 4.

#### **2. A Meta-Population Model on a Network**

The transmission of infectious diseases raises many important questions. In some instances, the average behavior of a large population with respect to the time is sufficient to provide useful insight from the available data. However, the spatial component of many transmission systems has been recognized to be of pivotal importance in the recent years. Due to this, spatially heterogeneous interventions must be included in the model, and hence it is essential to properly represent the transmission pattern. A reasonable hypothesis may consider that the spatial aspects of transmission heavily influence the aggregation characteristic of epidemic influence. Hence, we need to investigate data by using models that include such spatial connections. For example, the understanding of human mobility and the developing of qualitative and quantitative theories is of key importance for the modeling and for the comprehension of human infectious disease dynamics on geographical scales of different size.

#### *2.1. Spatial Heterogeneity in Epidemiological Models*

Ideally, the model should be able to account for the states of all *N* individuals in the population in an independent manner and, at the same time, it should allow for arbitrary interactions among them. The analysis of these models is a difficult task, and the computational cost of numerical simulations is very onerous and the extraction of the collective behaviors very complex. Although studies on the temporal dynamics of diseases proved insightful, incorporating space explicitly into epidemiological models revealed various emergent properties [19]. The phenomenon of the spatial spread of infection involves several components and scale [20]. Indeed, small region/group models can incorporate spatial heterogeneity, and more general models allowing for larger households with continuous or discrete time can be developed. Other typical approaches encompassing spatial variation in epidemic models involve partial differential equations (PDE). There exist, nonetheless, cases and scenarios where the latter type of spatial approach may not reliably model the phenomenon. Consider, for example, a human specific disease which is spread only by person-to-person contact and consider a geographical context consisting in a large country with a small number of large cities and a very sparse (or even non-existent) rural population [21]. The travel of individuals between discrete geographical regions and/or cities plays a pivotal role in the disease spreading. The depicted situation is easily described by a directed graph, where the vertices represent the cities (or discrete geographical regions/patches) and the arcs represent the links between such cities [22].

The main approaches for spatial models concern a different scale: an individual-based simulation, a meta-population model or a network model (see, e.g., [14]). Individual-based models explicitly represent every individual *i* with a state *Xi*(*t*) at time *t*, e.g., *Xi*(*t*) = 1 indicates that *i* is infected while *Xi*(*t*) = 0 indicates that *i* is healthy at time *t*. Infected nodes can transmit the disease to neighboring nodes and an individual becomes infected with a certain probability based on the status of neighboring individuals. The meta-population framework consists of dividing the whole population into distinct subpopulations, each having independent internal dynamics. In addition, at the same time, there is a limited interaction between the subpopulations. This approach has been used to great effect within ecological systems [23,24]. The individuals in such subpopulations belong to a particular state (e.g., susceptible or infectious) which can change during the time. For large networks, a general approach consists of merging some network information in a relatively small set of statistics and then studying the impact of such statistics on the infection spread [25].

In the following section, we describe our meta-population model, which considers communities as the aggregate unit that may represent either subpopulations in different areas or distinct groups with similar characteristics (e.g., students on a campus and citizens of a neighborhood or high school students or office workers). Then, each subpopulation is partitioned according to a particular state of individuals with respect to disease. Finally, connections and mobility between different communities are introduced. We point out that various proposed models encompass the geography of spread of the disease, but they do not present a mathematical analysis of their main properties, while presenting realistic simulations and an appropriate identification of the parameters involved, e.g., [26].

#### *2.2. A Prototype: SEIR Model on a Direct Graph*

We introduce a prototype model that can be generalized considering several states related to a given disease. Our analysis can, therefore, easily be extended to these more complex models. We partition a population of *N* individuals into subpopulations (groups, patches, communities, etc.) without taking into account any biological interpretation they have but considering spatially segregated large subpopulations. In this way, we can encompass a more realistic contact structure into epidemic models, since it usually preserves analytic tractability (in stochastic and in deterministic models), but at the same time it also captures the most important structural inhomogeneity in contact patterns in several applied contexts. The subpopulations and the interactions/connections between them are modeled through a weighted direct graph G = (V, E) with *n* vertices (nodes, regions, patches, subpopulations) and *m* edges (connections). Each edge is described by an ordered pairs of nodes (*u*, *v*), where *u*, *v* ∈ V. We label the nodes with an integer index; two vertices *i* and *j* of the directed graph are joined or adjacent if and only if there exists an edge from *i* to *j* or from *j* to *i*. If such an edge exists, then *i* and *j* are called its endpoints. If there is an edge from *i* to *j* then *i* and *j* are often called tail and head, respectively. The (*n* × *n*) *adjacency matrix Ad* associated to the graph is constructed as follows: if there exists an edge from node *i* to node *j*, then the entry at row *i* and column *j* is set to 1 in the matrix *Ad*: *ad ij* = 1.

In node *i*, the corresponding subpopulation possesses *Ni* individuals, and ∑*n i*=1 *Ni* = *N*. We hypothesize that individuals can move to a different node, interact with people in that node and then return to the original one. If *ad ij* = 1, there is an interaction between node *i* and node *j*, but not all the subpopulation *Ni* from node *i* interacts with the population in node *j*: we denote by *aij* the total amount of the subpopulation *i* that "goes" to node *j* and interacts with the people in that node. We call *A* the *routing matrix* with entries *aij*, so that ∑*n j*=1 *aij* = *Ni*, *i* = 1, ... , *n*. Associated to *A*, let *Po* the *probability outgoing matrix* with entries *po ij*, where we denote by *po ij* the percentage (probability) of the subpopulation *i* that "goes" to node *j*. In addition, we denote by *Pi* the *probability incoming matrix* with entries *pi ij*, where *pi ij* is now the percentage (probability) of the subpopulation in *j* that "arrived" from *i*. Finally, let *Mi* = ∑*n j*=1 *aji* be the total amount of people arrived in node *i* = 1, ... , *n*, so that ∑*n i*=1 *Mi* = *N* again. Then, for any *i* = 1, ... , *n*, ∑*n j*=1 *po ij* = ∑*n j*=1 *pi ji* = 1. Moreover we have

$$A = \text{Diag}(N\_1, N\_2, \dots, N\_n) \\ P^o = P^i \text{Diag}(M\_1, M\_2, \dots, M\_n)$$

where Diag(*x*1, *x*2, ... , *xn*) is the diagonal matrix with the vector (*x*1, *x*2, ... , *xn*)*T* ∈ R*n* on the main diagonal.

Four different discrete classes are considered for statuses of individuals in each: susceptible, exposed, infectious and recovered (SEIR model) [2]. All individuals are born as susceptible: a susceptible individual in contact with an infectious one may become exposed; the probability depends on the particular strain of the disease. Exposed individuals are infected but not yet infectious: individuals experience a long incubation duration. With a suitable incubation rate, latent individuals become infectious. Finally, a reliable assumption is that the immune system of infectious individuals combats the infection and then they move directly into the recovered class, which refers to individuals that are no longer infectious and have gained full immunity from further infection. Let *S*(*t*), *E*(*t*), *I*(*t*), *R*(*t*) the number of individuals in a node at time *t*, *S*(*t*) + *E*(*t*) + *I*(*t*) + *R*(*t*) = *N*: we consider a time interval in which we can neglect demographics. Without any interaction with other nodes, within a deterministic approach of the compartmental models, with continuous time *t*, the epidemic dynamics can be described by the system of differential equations in (1):

$$\begin{array}{l} \mathcal{S}(t) = -\lambda \, \mathcal{S}(t) \\ \dot{E}(t) = \lambda \, \mathcal{S}(t) - \mu E(t) \\ I(t) = \mu E(t) - \gamma I(t) \\ \mathcal{R}(t) = \gamma I(t) \end{array} \tag{1}$$

where the parameter *λ* is the force of infection, *γ* is the recovery rate and 1/*μ* is an average latent period.

With respect to the behavior of an epidemic, *λ* is the rate at which susceptible individuals become infected or exposed and it is a function depending on the number of infectious individuals; it contains information about the interactions between individuals that concur to the infection transmission. If we suppose that the population of *N* individuals mixes at random, meaning that all pairs of individuals have the same probability of interacting, the force of infection may be computed as:

$$\begin{array}{rcl} \lambda & = \text{transmission rate} \\ & \times \text{effective number of contacts per unit time} \\ & \times \text{proportion of contact infectious} \\ & \sim \tau \times n\_c \times \frac{I}{N} = \beta \frac{I}{N} \end{array}$$

Then the dynamics state,

$$\begin{aligned} \dot{S}(t) &= -\beta \frac{l}{N} S(t) \\ \dot{E}(t) &= \beta \frac{l}{N} S(t) - \mu E(t) \end{aligned} \tag{2}$$

where *β* is the infectious rate. Rescaling the quantities *S*, *E*, *I*, *R* dividing by *N* we obtain,

$$\begin{aligned} \dot{s}(t) &= -\beta \iota(t)s(t) \\ \dot{\iota}(t) &= \beta \iota(t)s(t) - \mu \iota(t) \\ \dot{\iota}(t) &= \mu \iota(t) - \gamma \iota(t) \\ \dot{\tau}(t) &= \gamma \iota(t) \end{aligned} \tag{3}$$

Pay attention to the fact that *ı*˙ stands for the derivative of the function *ı*. Now, we take a node *j* that is connected to the other nodes as encoded in matrix *A*. Then, *Sj*(*t*) can change due to the contribution of susceptible people from *j* that reached an adjacent node *k* and met infectious people in that node, wheresoever they came from. Then the contribution to *S*˙ *j* due to the interactions in node *k* is given by the *po jkSj* = *ajksj* susceptible people that met a population in node *k* with a proportion of infectious people given by

$$\begin{aligned} \frac{\#\{\text{inflection people in node } k\}}{\#\{\text{total people in node } k\}} &= \frac{\sum\_{l=1}^{n} p\_{lk}^{o} I\_{l}}{\sum\_{l=1}^{n} a\_{lk}} = \frac{\sum\_{l=1}^{n} p\_{lk}^{o} I\_{l}}{M\_{k}} \\ &= \frac{\sum\_{l=1}^{n} a\_{lk} I\_{l}}{\sum\_{l=1}^{n} a\_{lk}} = \sum\_{l=1}^{n} p\_{lk}^{\dot{i}} I\_{l}. \end{aligned}$$

Let the vectors **X**(*t*)=(*x*1(*t*), *x*2(*t*), ... , *xn*(*t*))*T* ∈ R*n*, with *X* = *S*, *E*, *I*, *R*, the *SEIR* model on the graph G is the following

$$\begin{aligned} \dot{S}(t) &= -\beta \text{Diag}(S(t)) \dot{\beta} I(t) & \dot{\mathbf{s}}(t) &= -\beta \text{Diag}(\mathbf{s}(t)) B \mathbf{r}(t) \\ \dot{E}(t) &= \beta \text{Diag}(S(t)) \dot{B} I(t) - \mu E(t) & \dot{\mathbf{e}}(t) &= \beta \text{Diag}(\mathbf{s}(t)) B \mathbf{r}(t) - \mu \mathbf{e}(t) \\ \dot{I}(t) &= \mu E(t) - \gamma I(t) & \dot{\mathbf{t}}(t) &= \mu \mathbf{e}(t) - \gamma \mathbf{r}(t) \\ \dot{R}(t) &= \gamma I(t) & \dot{\mathbf{r}}(t) &= \gamma \mathbf{r}(t) \end{aligned} \tag{4}$$

where *B*ˆ = *Po* Diag(*M*1, ... , *Mn*)−1*Po*, *B* = *PoPi*, and the equations on the right side have been obtained by a premultiplication with Diag(*N*1, *N*2,..., *Nn*)−1.

**Remark 1.** *We have assumed that the parameter β is the same in all nodes. It is possible to easily introduce a different parameter for each node considering more heterogeneity in the model.*

In the following, we adopt the notations **1** = [1, 1, ... , 1] , **0** = [0, 0, ... , 0] , for any vectors *x*, *y* ∈ R*n*, *x y* ⇔ *xi* < *yi*, *i* = 1, 2, ... , *n*; *x* ≤ *y* ⇔ *xi* ≤ *yi*, *i* = 1, 2, ... , *n* (and *x* < *y* if *x* ≤ *y* but *x* = *y*).

We suppose that the directed graph G is strongly connected, i.e., there exists a path in each direction between each pair of vertices of the graph, then the matrices *A* and *P* are irreducible. This means that we cannot divide the nodes of the graph into two subsets such that there are no connections between the nodes of the two subsets but only within each subset. It also follows that the matrix *B* is a non-negative irreducible *n* × *n* matrix; by the Perron–Frobenius theorem [27] we deduce:

Let *Bs*(*t*) = Diag(*s*(*t*)) *B*, *λs*(*t*) = *ρ*(*Bs*(*t*)) the dominant eigenvalue of *Bs*(*t*) and *vs*(*t*) 0 the corresponding positive left eigenvector.

It is easy to prove that this system of differential equations has a local solution by standard argument by the Cauchy–Lipschitz–Picard–Lindelöf theorem. Furthermore, if *ı*(0) > 0 and *s*(0) 0 then

Then, the solution is non-negative. Moreover, it is easy to check that if *sj*(0) + *ej*(0) + *ij*(0) + *rj*(0) = 1 then *sj*(*t*) + *ej*(*t*) + *ij*(*t*) + *rj*(*t*) = 1 and the solution is bounded, so there is a global solution for any time *t* > 0.

About the behavior of the epidemic dynamics, we will analyze the most important epidemiological properties:

**Theorem 1** (Asymptotic behavior)**.** *Denote by xs*(*t*) *the normalized version of vs*(*t*)*: xs*(*t*) *Bs*(*t*) = *λs*(*t*)*xs*(*t*) *and xs*(*t*)*xs*(*t*) = 1*. For any initial condition, λs*(*t*) *is a continuous function, and there exist the limits of the above quantities: s*∞ = lim*t*→∞ *s*(*t*)*,* lim*t*→∞ *Bs*(*t*) = *Bs*∞ = Diag(*s*∞)*B,* lim*t*→∞ *λs*(*t*) = *λs*∞ *, where ρ*(*Bs*∞ ) = *λs*∞ *.*

*If s*∞ > 0*, then ρ*(*Bs*∞ ) = *λ*∞ > 0 *and any converging subsequence of xs*(*t*) *converges to a λ*∞*-eigenvector.*

*If, in addition, s*∞ 0*, λ*∞ *is simple, and then* lim*t*→∞ *xs*(*t*) = *xs*∞ 0*, where x s*∞ *Bs*∞ = *λs*∞ *x s*∞ *.*

**Proof.** The existence of the limit of *s*(*t*) is obvious, since *s*(*t*) is a continuous monotone function. Accordingly, *Bs*(*t*) = Diag(*s*∞)*B* converges to Diag(*s*∞)*B*. From now on, to simplify the notations in the proof, we will use *Bt* = *Bs*(*t*), *λt* = *λs*(*t*) and *xt* = *xs*(*t*) for any *t* ∈ R+. Now, let *t* fixed and *δ* > 0 sufficiently small so that (1 − ) Diag(**1**) ≤ Diag(*sv*) Diag(*st*)−1 ≤ (1 + ) Diag(**1**) for *v* ∈ (*t* − *δ*, *t* + *δ*). Then, for *t* ≤ *v* < *t* + *δ*,

$$\begin{split} 0 \le \lambda\_t - \lambda\_v \le \lambda\_t - \min\_i \frac{[\mathbf{x}\_t^\top B\_v]\_i}{[\mathbf{x}\_t]\_i} \\ &= \lambda\_t - \min\_i \frac{[\mathbf{x}\_t^\top \text{Diag}(\mathbf{s}\_v) B]\_i}{[\mathbf{x}\_t]\_i} \\ &= \lambda\_t - \min\_i \frac{[\mathbf{x}\_t^\top \text{Diag}(\mathbf{s}\_v) \text{Diag}(\mathbf{s}\_t)^{-1} \text{Diag}(\mathbf{s}\_t) B]\_i}{[\mathbf{x}\_t]\_i} \\ &\le \lambda\_{t1} - \min\_i (1 - \epsilon) \frac{[\mathbf{x}\_t^\top (B\_t)]\_i}{[\mathbf{x}\_t]\_i} \\ &\le \epsilon \lambda\_{t\_t} \end{split}$$

while, for *t* − *δ* < *v* ≤ *t*,

$$\begin{aligned} 0 \le \lambda\_v - \lambda\_t &\le \max\_i \frac{[\mathbf{x}\_t^\top B\_v]\_i}{[\mathbf{x}\_t]\_i} - \lambda\_t \\ &= \max\_i \frac{[\mathbf{x}\_t^\top \text{Diag}(\mathbf{s}\_v) \text{Diag}(\mathbf{s}\_t)^{-1} B\_t]\_i}{[\mathbf{x}\_t]\_i} - \lambda\_t \\ &= \max\_i (1 + \epsilon) \frac{[\mathbf{x}\_t^\top (B\_t)]\_i}{[\mathbf{x}\_t]\_i} - \lambda\_t \\ &\le \epsilon \lambda\_{t\prime} \end{aligned}$$

which implies that *λs*(*t*) is a continuous monotone function that must have a limit; denote it by *λ*∞. If *s*∞ = **0**, then **1***Bt* → **0**, which implies that *λt* → 0. From now on, we then assume **1***s*∞ = |*s*∞|1 > 0, so that *λ*∞ > 0. Let *xtn* any converging subsequence, call *x*˜ its limit. Then, is a *λ*∞ non-negative eigenvector of *B*∞, since

$$\begin{split} \boldsymbol{\tilde{\mathfrak{x}}}^{\top} \boldsymbol{B}\_{\infty} &= \bullet \big( \boldsymbol{\tilde{\mathfrak{x}}}^{\top} - \boldsymbol{\mathfrak{x}}\_{t\_{n}} \big) \boldsymbol{B}\_{\infty} + \boldsymbol{\mathfrak{x}}\_{t\_{n}}^{\top} \big( \boldsymbol{B}\_{\infty} - \boldsymbol{B}\_{t\_{n}} \big) + \boldsymbol{\mathfrak{x}}\_{t\_{n}}^{\top} \boldsymbol{B}\_{t\_{n}} \\ &= \boldsymbol{\varepsilon}(\boldsymbol{n}) \boldsymbol{\mathfrak{1}}^{\top} + \boldsymbol{\lambda}\_{t\_{n}} \boldsymbol{\mathfrak{x}}\_{t\_{n}}^{\top} \boldsymbol{B}\_{t\_{n}} \\ &\to \boldsymbol{\lambda}\_{\infty} \boldsymbol{\mathfrak{x}}^{\top} .\end{split}$$

Now, if we add the hypothesis that *s*∞ 0, then *B*∞ is still a Perron matrix, and hence there exists a unique positive eigenvector of *B*∞, whence *xt* → *x*∞ 0.

**Theorem 2** (Threshold)**.** *Consider the SEIR model* (4) *on a strongly connected graph G, let λs*(*t*) *be the dominant eigenvalue of Bs*(*t*) *and let vs*(*t*) 0 *be the corresponding positive left eigenvector.*

**Proof.** Take ∈ 0, *γ*−*βλs* (*τ*) *βλs* (*τ*) , so that

$$
\beta \lambda\_{\mathfrak{s}}(\tau)(1+\epsilon) - \gamma < 0,
$$

and define

$$c\_{\mathfrak{e}} = \min\left(\gamma - \beta \lambda\_{\mathfrak{s}}(\tau)(1+\mathfrak{e}), \frac{\mathfrak{e}\mu}{1+\mathfrak{e}}\right) > 0. \tag{5}$$

Multiplying the weighted sum of the second equation of exposed and the third equation of the infected by *vs*(*τ*),

$$\begin{split} \mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})^{\top}(\dot{\mathfrak{e}}(t)(1+\mathfrak{e})+\mathfrak{i}(t)) &= \mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})^{\top}(\beta \operatorname{Diag}(\mathfrak{s}(t))\operatorname{Br}(t)(1+\mathfrak{e}) - \mathfrak{e}\mu\mathfrak{e}(t) - \gamma\mathfrak{i}(t)) \\ &= \mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})^{\top}(\beta B\_{\mathfrak{s}}(t)\mathfrak{i}(t)(1+\mathfrak{e}) - \mathfrak{e}\mu\mathfrak{e}(t) - \gamma\mathfrak{i}(t)), \end{split}$$

then, for *t* ≥ *τ*,

$$\frac{d\mathfrak{v}\_{\mathfrak{s}}(\pi)^{\top}(\mathfrak{e}(t)(1+\mathfrak{e})+\mathfrak{r}(t))}{dt} \leq \mathfrak{v}\_{\mathfrak{s}}(\pi)^{\top}(\beta B\_{\mathfrak{s}}(t)(1+\mathfrak{e})\mathfrak{r}(t) - \mathfrak{e}\mu\mathfrak{e}(t) - \gamma\mathfrak{r}(t)).$$

$$\text{where } \mathfrak{e} \text{ and } \mathfrak{r} \text{ are } \omega\text{-a) } \mathfrak{T}\_{\mathfrak{s}}\mathfrak{a} \text{ and } \mathfrak{r} \text{ are }$$

Now, *βvs*(*τ*)*Bs*(*τ*) = *vs*(*τ*)*βλs*(*τ*), then

$$\begin{aligned} \frac{d\,\mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})\,^{\top}(\mathfrak{e}(t)(1+\mathfrak{e})+\mathfrak{r}(t))}{dt} & \leq \mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})^{\top}((\mathfrak{\beta}\lambda\_{\mathfrak{s}}(\mathfrak{r})(1+\mathfrak{e})-\gamma)\mathfrak{r}(t)-\mathfrak{e}\mu\mathfrak{e}(t)) \\ & \leq -\mathfrak{v}\_{\mathfrak{s}}(\mathfrak{r})^{\top}\mathfrak{c}\_{\mathfrak{e}}(\mathfrak{r}(t)+(1+\mathfrak{e})\mathfrak{e}(t)), \end{aligned}$$

where *c* is defined in (5). Using the previous differential inequality, Gronwall lemma implies that

$$
\mathfrak{w}\_{\mathfrak{s}}(\mathfrak{r})^{\top}(\mathfrak{e}(t)(1+\mathfrak{e})+\mathfrak{r}(t)) \leq \mathfrak{w}\_{\mathfrak{s}}(\mathfrak{r})^{\top}(\mathfrak{e}(\mathfrak{r})(1+\mathfrak{e})+\mathfrak{r}(\mathfrak{r}))e^{-\mathfrak{c}\_{\mathfrak{e}}(t-\mathfrak{r})}.
$$

To prove (2), we start by noticing that

$$\frac{d}{dt} \left( \mathfrak{v}\_{\mathfrak{s}}(0)^\top \left( \mathfrak{e}(t) + \mathfrak{e}(t) \right) \right)\_{t=0} = \left( \mathfrak{f} \lambda\_{\mathfrak{s}}(0) - \gamma \right) \mathfrak{v}\_{\mathfrak{s}}(0)^\top \mathfrak{e}(0).$$

Since (*βλs*(0) − *γ*) > 0, *ı*(0) > **0**, and *vs*(0) **0**, then

$$\frac{d}{dt}\left(\mathfrak{v}\_{\mathfrak{s}}(0)^{\top}(\mathfrak{e}(t)+\mathfrak{r}(t))\right)\_{t=0} > 0.$$

We have that the solution is a continuous differentiable function, then exists *t* ∗ > 0 such that for *s* ∈ (0, *t* ∗)

$$\frac{d}{dt} \left( \mathfrak{v}\_{\mathfrak{s}} \left( \mathfrak{v} \left( \mathfrak{d} \right)^{\top} \left( \mathfrak{e} \left( t \right) + \mathfrak{e} \left( t \right) \right) \right)\_{t=\mathsf{s}} > 0, \mathsf{d} \right)$$

which implies that *q*0(*t*) = *vs*(0) (*e*(*t*) + *ı*(*t*)) increases for *t* ∈ (0, *t* ∗).

Moreover, it is *s*˙ **0** and *s*(*t*) **0**, then exist lim*t*→+∞ *s*(*t*) = *s*∞. Define *b j* the *j*–th row of *B*, so that *s*˙*j* = *βsjb j ı*. Since

$$\ddot{s}\_{\dot{j}} = \beta^2 s\_{\dot{j}} (\boldsymbol{b}\_{\dot{j}}^\top \boldsymbol{\nu})^2 + \beta s\_{\dot{j}} \boldsymbol{b}\_{\dot{j}}^\top (\boldsymbol{\mu}\boldsymbol{e} - \boldsymbol{\gamma}\boldsymbol{\nu})$$

then *s*¨ ≤ *K*, and hence the boundedness and monotonicity of *s*(*t*) implies lim*t*→+∞ *s*˙(*t*) = **0**, which implies **1***s*˙(*t*) → 0. Then, for any > 0 we have for *t* sufficiently large that

$$\begin{aligned} \frac{d}{dt} \mathbf{1}^\top \mathbf{e}(t) &\leq \frac{d}{dt} \Big( \mathbf{1}^\top (\mathbf{e}(t) + \mathbf{s}(t)) \Big) + \epsilon \\ &\leq \epsilon - \mu \mathbf{1}^\top \mathbf{e}(t). \end{aligned}$$

Then *e*(*t*) → **0**, that also implies **1***e*˙(*t*) → 0. As a consequence, for *t* sufficiently large,

$$\begin{aligned} \frac{d}{dt} \mathbf{1}^\top \mathfrak{a}(t) &\leq \frac{d}{dt} \Big( \mathbf{1}^\top (\mathfrak{e}(t) + \mathfrak{s}(t) + \mathfrak{e}(t)) \Big) + \mathfrak{e} \\ &\leq \mathfrak{e} - \gamma \mathbf{1}^\top \mathfrak{a}(t) \end{aligned}$$

so that *ı*(*t*) → **0**.

#### *2.3. A Control Problem*

We have adopted a network framework that explicitly accounts for the interactions structure among individuals and group of individuals, in order to provide insights regarding the spread of a disease. If the proposed model describes the epidemiological phenomenon sufficiently well, some problems relating to the behavior and the forecast of the epidemic itself can be addressed.

First, we would like to prevent an epidemic. This is achieved when condition (1) in Theorem 2 holds at *t* = 0. Before the epidemic starts, the fractions of infected/exposed individuals are negligible, for viral infections the recovery rate *γ* is usually out of control. Then, the only way to satisfy the no-epidemic requirement is either: control the transmission (which means to reduce *β* and/or interactions) or immunization (meaning to increase *r*(0)). Second, we aim to limit the economic and social impact as the epidemic occurs. The supply of healthcare services is inelastic in the short run. Thus, it is important to maintain the maximum infection rate below the capacity of the existing healthcare system. This may be achieved by lowering the transmission rate, by controlling the inflow and the outflow of individuals from and into a node.

We point out that only recent works, e.g., [28–30], started investing the trade-off between epidemic and economic costs with some analysis. The aim of our applications would like to be a new step in this direction inside a well based framework. We introduce the following diagonal matrix

$$\mathcal{U} = \text{Diag}(\mu\_{\text{loc}(i)})\_{i=1\prime}^{n} \quad \mathcal{V} = \text{Diag}(\upsilon\_{\text{loc}(i)})\_{i=1\prime}^{n} \tag{6}$$

where *uloc*(*i*) ∈ (0, 1], *i* = 1, ... , *n* are the control variables for the incoming individuals into the node *i*, while *vloc*(*i*) ∈ [0, 1], *i* = 1, ... , *n* are the control variables of the outgoing individuals from the node *i* to other nodes. Then, the routing matrix *A*, and its associated matrices (see their definitions in Section 2.2), changes as follows

$$\begin{aligned} \bar{A}\_{uv} &= \mathsf{U}A \, V = \text{Diag}(\hat{\mathsf{N}}\_{1}, \dots, \hat{\mathsf{N}}\_{n}) \hat{P}^{\hat{o}} = \bar{P}^{\hat{i}} \, \text{Diag}(\hat{\mathsf{M}}\_{1}, \dots, \hat{\mathsf{M}}\_{n}), \\ \hat{\tilde{B}} &= \bar{P}^{\hat{o}} \, \text{Diag}(\tilde{\mathcal{M}}\_{1}, \dots, \tilde{\mathcal{M}}\_{n})^{-1} \bar{P}^{\hat{o}^{\top}}, \qquad \bar{B} = \bar{P}^{\hat{o}} \bar{P}^{\hat{i}^{\top}}, \end{aligned}$$

so that the *SEIR* model on the "controlled" graph G becomes

$$\begin{aligned} \dot{S}(t) &= -\beta \operatorname{Diag}(\mathcal{S}(t)) \tilde{\beta} I(t) & \dot{\mathfrak{s}}(t) &= -\beta \bar{I} \operatorname{Diag}(\mathfrak{s}(t)) B \mathfrak{r}(t) \\ \dot{E}(t) &= \beta \operatorname{Diag}(\mathcal{S}(t)) \tilde{\beta} I(t) - \mu \mathcal{E}(t) & \dot{\mathfrak{s}}(t) &= \beta \bar{I} \operatorname{Diag}(\mathfrak{s}(t)) B \mathfrak{r}(t) - \mu \mathfrak{e}(t) \\ \dot{I}(t) &= \mu \mathcal{E}(t) - \gamma I(t) & \dot{\mathfrak{r}}(t) &= \mu \mathfrak{e}(t) - \gamma \mathfrak{r}(t) \\ \mathcal{R}(t) &= \gamma I(t) & \dot{\mathfrak{r}}(t) &= \gamma \mathfrak{r}(t) \end{aligned}$$

where ˜*I* = Diag(*N*˜ 1/*N*1, ..., *N*˜ *n*/*Nn*).

In order to study the "lockdown policies" applied to the various groups (nodes), we combine a measure of social cost (i.e., hospitalization cost) and, above all, loss of life and the economic loss. The first objective consists of minimizing both total (excess) deaths during the epidemic and the public health cost. We suppose that it is possible estimating the severity of the epidemic in a time interval [0, *T*], by weighing the total number of infected

$$\text{social cost} = \text{C}^T \int\_0^T I(t)dt$$

where *C* is a vector of positive weights, and the integral is to be understood component by component. The lockdown of individuals affects the economic activities; we model the economic loss through the evaluation of the reduction in flow of individuals between nodes with a linear cost function

$$\text{economic loss} = \mathbf{W}^T \int\_0^T \left[ A - \tilde{A}\_{\mu\upsilon} \right] \cdot \mathbf{1} dt.$$

The goal is finding an optimal trade-off between the total economic loss and the total social cost. Then, the optimal strategy is obtained by minimizing the following cost function

$$\min\_{\mathbf{U},\mathbf{V}} \left( \mathbf{C}^T \int\_0^T \mathbf{I}(t)dt - \mathbf{W}^T \int\_0^T \left[A - \tilde{A}\_{uv}\right] \cdot \mathbf{1} dt \right), \tag{7}$$

where *T* is the time for which a certain strategy is applied.

**Remark 2.** *If the meta-population model represents a non-geographical subdivision, but instead it is dependent on certain individuals' characteristics (such as age, profession, habits, etc.), the weights C and W can include information linked to these characteristics (e.g., propensity for mobility, disease mortality). In this case, the lockdown strategy can change based on the vulnerability of each group.*

#### **3. Numerical Tests**

This section is devoted to numerical experiments. These experiments were carried out on a laptop equipped with Linux 19.04, with an Intel(R) Core(TM) i5–8250U CPU (1.60 GHz), 16 GiB RAM memory (Intel, Santa Clara, CA, USA) and under MATLAB R2020b environment (MathWorks, Natick, MA, USA).

In order to use our framework, two types of parameters are needed:

For the first set of parameters, we have referred to a recent work [31] in which the authors applied a SEIR epidemiological model to the recent SARS-CoV-2 outbreak in the world. Moreover, they focused on the application of a stochastic approach in fitting the biological model parameters analyzing the official data and the predicted evolution of the epidemic in the Italian regions, Spain and South Korea. We considered two different scenarios,

For the topology of the directed graph G(*V*, *E*), we did not refer to any particular geographic area but we reproduced a realistic situation. The network consists of three large agglomerates, each one representing a city. The nodes of the graph are the neighborhoods of the cities and the edges represent the connections between such neighborhoods: these edges encompass the social and working movements between the nodes; hence, they are not simply geographical connections (see Section 2.2). The number of the nodes is 20, 10 and 5, respectively, meaning that the largest city has 20 neighborhoods, the second

one 10 and the last one just 5 neighborhoods. This toy model considers a social cost *C* which is ten times higher than the economic cost *W*: normalizing such costs leads to set *C* = **1**,*W* = 10 · **1**. The matrix *A* is set starting from the adjacency matrix *E* and from the population of the nodes: the number of individuals, the subpopulations and the matrices *Po* and *Pi* were randomly selected using suitable probability distributions.

All scenarios started with the initial distribution for susceptible, infected and exposed individuals: the epidemic starts from 1/5 of randomly selected neighborhoods of the largest city. Once a lockdown strategy is decided by optimizing Equation (7), it is applied for 14 days: after such period the new distributions for *S*, *I* and *E* are checked and a new optimization is carried on. The last time interval has a longer duration for observing the effects of the overall strategy on the long period. This scheme is applied three times in the numerical simulations.

Figure 1a presents the optimized strategy for Scenario A, while Figure 1b shows the strategy for Scenario B. We have not represented the whole network but a part of it considering nodes that represent agglomerations with a different number of individuals. Furthermore, the strategy that optimizes our objective function is reproduced by showing the values of the vector *V*, see (6), for some cities/agglomerations present in the network. Values close to 1 mean that there are no particular restrictions on mobility, values close to 0 mean strong restrictions on movement. We point out that the value 0 is not allowed because it is not realistic to consider a total block of each movement in this context.

In the case of Scenario A, the parameter of the disease induced a light lockdown (85%) on the large city (blue line), whilst the other two are almost completely open. On the second case, the disease is more infectious: the large city is forced to adopt a severe lockdown, while the strategy on the other two suggests a mild lockdown. As soon as the epidemic spreads, due to the characteristics of the disease, even the smallest cities are forced to adopt a more severe strategy until the number of infected individuals decreases. We can observe in Scenario B that, after a period of severe lockdown, when the last chosen strategy is applied for a longer period, then a further severe approach must be adopted in order to contain the epidemic.

In both scenarios, we can observe that there is a converging behavior of the strategies to be applied on the three different cities. Check, for example, the period 50–100 days in Scenario A: the lockdown strategy for the largest city (blue line) is increasing, while the strategies for the other cities (orange and yellow lines) are decreasing. Eventually, due to the disease parameters, they converge to 1, meaning that the epidemic threat is no more. This behavior is more evident in Scenario B.

**Figure 1.** (**a**,**b**): Lockdown strategies for Scenario A and B, respectively. We show the optimal strategy *V* that minimizes (7) the closer is to 0, the more severe the restrictions are. On the other hand, values close to 1 denote very mild restrictions on mobility.

Figure 2 presents a visual representation of the evolution of the strategies for Scenario B at *t* = 0, 20, 40 and 179 days. The color of the nodes represents the number of infected people in that node, while the color of the edges represents the percentage of people blocked. At the beginning, the optimal strategy suggests to mainly block the outgoing from the large city (the agglomerate on the top left corner), while the connection between and inside the other two are open. Letting the epidemic spread in a controlled way, in order to maintain economy, induces an increase in the number of infected people, then the connections between cities must be reduced (the more red the edges are, the less people are allowed to move).

**Figure 2.** (**a**–**d**): Evolution of the strategy and infected people in Scenario B for timestamps *t* = 0, 20, 40, 179, respectively. The color of the nodes represents the number of infected people in that node, while the color of the edges represents the percentage of people blocked.

#### **4. Conclusions**

The network perspective allows to relax the assumption of uniform random mixing and then we are able to model the population interaction patterns during epidemics. Moreover, a network-based model provides useful and important insights about the spread of a disease; such insights cannot be inferred using the classical model. We have focused on a SEIR meta-population model on a network in order to characterize the epidemic dynamics and to predict possible contagion scenarios.

A network model can encompass differences in the number of interactions between individuals in a population. For example, there exist some realities where people may live in small environments and they have relatively few contacts (work and/or social life), while others may live in dense and more populated centers, where the usage of public and maybe crowded public transportation is very common, they work in high-contact environments and they have a large number of interactions with many others outside of work. The classical SEIR model does not allow to include such heterogeneity, while a network model can easily encompass it. Furthermore, it is possible to adapt our model and, instead of a geographical distinction of the subgroups of individuals, different stratifications of the population can be considered.

We justified the model by introducing and analyzing some of its properties, in particular, we proved a threshold theorem involving both biological parameters and the topology of the network. In a future paper, we will consider both time-dependent parameters and a detailed analysis of the asymptotic behavior of the solution of the proposed model. We point out that our analysis can be applied to recent models, e.g., [26], where numerical and statistical but not analytical results are provided.

Only recent works, e.g., [28–30], started investigating the trade-off between epidemic and economic costs with some analysis. In order to take a new step in this direction, we have also identified an optimal control problem that considers the advantages and benefits that arise from the application of optimal targeted policies, which lock down the various groups in an inhomogeneous way. The focus is on the balance between economic loss and loss of life. The main economic damages consist of lost productivity due to illness and also in the forgone productivity contributions of the blocked subpopulations. The lives lost are estimated via the number of the infected individuals, assuming that these losses represent a constant percentage of the latter.

Some preliminary numerical tests are provided; more effective results could be obtained by considering a suitable fitting of the parameters and based on some particular topology of the network. These issues will be analyzed in a future paper involving a different source of data [32] and recent optimization tools [33,34].

**Author Contributions:** Conceptualization, G.A. and G.N.; methodology, G.A., A.B. and G.N.; software, G.A.; validation, G.A. and A.B.; formal analysis, G.A., A.B. and G.N.; investigation, G.A. and G.N.; resources, G.A., A.B. and G.N.; data curation, G.A. and A.B.; writing—original draft preparation, G.A., A.B. and G.N.; writing—review and editing, A.B.; visualization, A.B.; supervision, G.A., A.B. and G.N.; project administration, G.N.; funding acquisition, G.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding from PRECISION project.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Acknowledgments:** We acknowledge people from ADAMSS center for the insights on the developing of the model.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**