**1. Introduction**

In the early 2020s, the spread of the COVID-19 pandemic highlighted the close link between economics and health in the context of emergency management. Because of this, assessing the impact of an epidemic phenomenon on a country's economy has emerged as one of the key aspects to consider in the context of containment strategies. From a mathematical point of view, a systematic approach to the study of the effects on the economies of countries facing a severe pandemic is a very complex problem and a mathematical model can only provide rough indications of the possible consequences, based on simplifying assumptions about the key parameters driving the pandemic evolution. The basic idea is to trace these phenomena back to the evolution of the so-called wealth distribution of a country, which measures how many people belong to increasing income levels.

A first attempt to understand changes in wealth distribution in the presence of epidemic spread was proposed in [1] by combining the classical SIR compartmental model of susceptible, infected and recovered individuals [2,3] with the kinetic model of wealth distribution introduced in [4], and assuming that, due to the presence of the pandemic, individuals in different compartments act differently in the economic process. Although the model was developed in a relatively simplified context, it has provided a general framework for socio-epidemiological modeling that can be easily extended to more complex dynamics, both in terms of economic transactions [5] and in terms of epidemic interactions [6,7]. We

**Citation:** Bernardi, E.; Pareschi, L.; Toscani, G.; Zanella, M. Effects of Vaccination Efficacy on Wealth Distribution in Kinetic Epidemic Models. *Entropy* **2022**, *24*, 216. https://doi.org/10.3390/e24020216

Academic Editors: Ryszard Kutner, H. Eugene Stanley and Christophe Schinckus

Received: 7 January 2022 Accepted: 25 January 2022 Published: 29 January 2022

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

**Copyright:** © 2022 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/).

mention in this direction the recent survey reported in [8] and the seminal approaches proposed in [9–12] investigating the economic effects of infectious diseases, as well as the study presented in [13].

More precisely, according to [4], the financial transactions in [1] were based on the choice of two parameters. The first defines the so-called safeguard threshold, i.e., the maximum percentage of money that the individual is willing to employ in a transaction, and the second is the random risk inherent in the transaction, characterized by its variance through a spread proportional to the square of the individual's wealth. There, the time dependence of the variance was postulated by assuming that, in the presence of a significant epidemic spread, the variance of the risk tends to increase. This is in agreemen<sup>t</sup> with the financial market reactions that were often observed during the COVID-19 pandemic to announcements of rising numbers of infected people in several countries [14]. With the use of the model in [1], it was possible to qualitatively observe the effects of the pandemic in terms of a reduction in the middle class and the increase in social inequalities (see also [15,16]).

The possibility, starting in early 2021, of launching a widespread vaccination campaign has led to general optimism about the ability to improve economic performance while limiting the health consequences of the epidemic. However, it is clear that the reduction of economic consequences is closely linked to the effectiveness of the vaccine in containing infections.

In this paper we will focus, at the level of wealth distribution, on the economic improvements induced by the vaccination campaign in terms of its percentage of effectiveness. The interplay between the economic trend and the pandemic will be evaluated by resorting to a mathematical model joining a kinetic model of wealth distribution based on binary transactions with a compartmental epidemic model including vaccinated individuals (see also [17]). In particular, a fraction of vaccinated individuals, which is determined by the efficacy of the vaccine, may contract the disease. Without intending to review the extensive literature on this topic, we cite the recent papers [18–26] that highlight the possible partial immunity provided by vaccinations. Moreover, the emergence of viral variants means that the efficacy of the vaccine inherently non-constant and subject to collective compliance with non-pharmaceutical interventions.

The underlying theoretical framework we consider is that of kinetic models for collective social phenomena, which allows for the linking of microscopic agent-based behavior to emerging observable patterns [27]. In particular, mathematical modeling of wealth distribution has seen a marked development in recent decades [5,28–35], in which, at least partially, the essential economic mechanisms that are responsible for the formation of large-scale economic indicators such as the Pareto or Gini index have been understood [36,37].

The interplay between epidemic spread and the social economic background is described here as the result of interactions among a large number of individuals, each of which is characterized by the variable *w* ∈ R+, measuring the amount of wealth of a single agent. In this regard, as shown in [1,8,38,39], the fundamental tools of statistical physics allow the understanding of epidemiological dynamics by linking classical compartmental approaches with a statistical description of economic aspects. Indeed, the multiscale nature of kinetic theory allows for the determination of the macroscopic (or aggregate) and measurable features of disease evolution [27,40,41].

The rest of the paper is organized as follows. Section 2 introduces the SIR-type system of kinetic equations that includes vaccinated individuals and combines the dynamics of wealth evolution with the spread of infectious disease in a system of interacting agents. Next, in Section 3 we study the main mathematical properties of the system, and show that, through a suitable asymptotic procedure, the solution of the kinetic system tends to the solution of a system of Fokker–Planck-type equations, which exhibits explicit equilibria of the inverse Gamma type. Finally, in Section 4, we investigate numerically the solutions of the Boltzmann-type kinetic system, and its Fokker–Planck asymptotics, along with the evolution of the Gini index, characterizing the wealth inequalities. These simulations

confirm the model's ability to describe phenomena that are characteristic of economic trends in situations compromised by the rapid spread of an epidemic, and their variations as a function of the effectiveness of the vaccination campaign.

#### **2. Wealth Dynamics in Epidemic Phenomena**

In this Section, we present an extension of the SIR-kinetic compartmental description of epidemic spreading introduced in [1], which additionally takes into account the population of vaccinated individuals. The model consists of a system of four kinetic equations describing the evolution of wealth in the presence of an infectious disease with partial efficacy of vaccination. The entire population is divided then into four compartments: susceptible individuals (*S*), who can contract the disease; identified infectious individuals (*I*), who are recognized to have contracted the disease and can transmit it; vaccinated individuals (*V*), who have received a vaccine, but can still be at least partially infected and contagious; and the recovered individuals (*R*), who are healed and immune. The model can be easily adapted to include disease-related mortality and other compartments of interest in terms of available data, such as records of hospitalized individuals. We refer to [3,6,7,42] and the references therein for possible developments in these directions. It should be noted that, since we are referring to an advanced epidemic situation in which we assume the existence of a vaccine, the dynamics of unidentified asymptomatic individuals, so significant in the early stages of the COVID-19 pandemic, has become less relevant thanks to mass screening programs. For this reason, we have chosen to employ only one compartment *I* related to the identified infected individuals. To measure the aggregate effects of vaccination over the whole population, we have considered the compartment *V* with a given vaccine efficacy.

The agents of each compartment are characterized uniquely by their wealth *w* ≥ 0. Hence, we denote by *fH*(*<sup>w</sup>*, *t*), *H* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}, the distributions of wealth at time *t* ≥ 0 in each compartment, such that *fH*(*<sup>w</sup>*, *t*)*dw* denotes the fraction of agents belonging to the compartment *J*, which, at time *t* ≥ 0, are characterized by wealth between *w* and *w* + *dw*. The total wealth distribution density is then defined by the sum of the distributions in all compartments

$$f(w,t) = f\_{\mathcal{S}}(w,t) + f\_{\mathcal{I}}(w,t) + f\_{\mathcal{V}}(w,t) + f\_{\mathcal{R}}(w,t), \qquad \int\_{\mathbb{R}\_+} f(w,t)dw = 1, \ \mathcal{S}$$

for all *t* ≥ 0. Hence, the fractions of the population belonging to each compartment are given by

$$J(t) = \int\_{\mathbb{R}\_+} f\_I(w\_\prime t) dw\_\prime \qquad J \in \{\mathcal{S}, I, V, \mathcal{R}\}.$$

We denote by *mJ*,*<sup>κ</sup>*(*t*) the local momenta of order *κ* for the wealth distributions in each compartment

$$m\_{\kappa,l}(t) = \frac{1}{f(t)} \int\_{\mathbb{R}\_+} w^\kappa f\_l(w, t) dw,\tag{1}$$

and we denote with *<sup>m</sup>κ*(*t*) the moment of order *κ* > 0 of the wealth distribution *f*(*<sup>w</sup>*, *t*)

$$m\_{\mathbb{k}}(t) = \int\_{\mathbb{R}\_+} w^{\mathbb{k}} f(w, t) dw = \sum\_{J \in \{\mathbb{S}, I, V, \mathbb{R}\}} J(t) \, m\_{\mathbb{k}, I}(t).$$

#### *2.1. The Kinetic Model*

Following [1], we assume that the evolution of the densities obeys an SIR-type compartmental model and that the wealth exchange process is influenced by the epidemic's dynamics. This gives a system of four kinetic equations for the unknown distributions *fH*(*<sup>w</sup>*, *t*), *H* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}, expressed by

$$\begin{aligned} \partial\_{l} f\_{\mathcal{S}}(\boldsymbol{w},t) &= -K(f\_{\mathcal{S}},f\_{l})(\boldsymbol{w},t) - a f\_{\mathcal{S}}(\boldsymbol{w},t) + \sum\_{\boldsymbol{I} \in \{\mathcal{S},\boldsymbol{I},\boldsymbol{V},\boldsymbol{R}\}} Q\_{\mathcal{S}}(f\_{\mathcal{S}},f\_{\boldsymbol{I}})(\boldsymbol{w},t), \\ \partial\_{l} f\_{\mathcal{I}}(\boldsymbol{w},t) &= K(f\_{\mathcal{S}},f\_{\mathcal{I}})(\boldsymbol{w},t) + (1-\zeta)K(f\_{\mathcal{V}},f\_{\mathcal{I}})(\boldsymbol{w},t) - \gamma\_{\mathcal{I}}f\_{\mathcal{I}}(\boldsymbol{w},t) + \sum\_{\boldsymbol{I} \in \{\mathcal{S},\boldsymbol{I},\boldsymbol{V},\boldsymbol{R}\}} Q\_{\mathcal{I}}(f\_{\mathcal{I}},f\_{\boldsymbol{I}})(\boldsymbol{w},t), \\ \partial\_{l} f\_{\mathcal{V}}(\boldsymbol{w},t) &= a f\_{\mathcal{S}}(\boldsymbol{w},t) - (1-\zeta)K(f\_{\mathcal{V}},f\_{\mathcal{I}})(\boldsymbol{w},t) + \sum\_{\boldsymbol{I} \in \{\mathcal{S},\boldsymbol{I},\boldsymbol{V},\boldsymbol{R}\}} Q\_{\mathcal{V}}(f\_{\mathcal{V}},f\_{\boldsymbol{I}})(\boldsymbol{w},t), \\ \partial\_{l} f\_{\mathcal{R}}(\boldsymbol{w},t) &= \gamma\_{\mathcal{I}}f\_{\mathcal{I}}(\boldsymbol{w},t) + \sum\_{\boldsymbol{I} \in \{\mathcal{S},\boldsymbol{I},\boldsymbol{V},\boldsymbol{R}\}} Q\_{\mathcal{R}}(f\_{\mathcal{R}},f\_{\boldsymbol{I}})(\boldsymbol{w},t), \end{aligned} \tag{2}$$

where *γ* ≥ 0 is the recovery rate for the infected compartment and *α* ∈ [0, 1] is the vaccination rate of individuals, whereas the term 0 ≤ 1 − *ζ* ≤ 1 quantifies the effectiveness of the vaccine, in such a way that high effectiveness corresponds to values close to one of the parameters, *ζ*. The operator *<sup>K</sup>*(·, ·) governs the transmission of the infection and is considered to be of the following form

$$\mathcal{K}(f\_{H\ast}f\_I)(w,t) = f\_H(w,t) \int\_{\mathbb{R}\_+} \beta(w,w\_\ast) f\_I(w\_\ast,t) \, dw\_\ast \, . \tag{3}$$

for any *H* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}. In (3) the function *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) ≥ 0 denotes the contact rate between people with wealth *w* and, respectively, *w*<sup>∗</sup>. A leading example for *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) is obtained by choosing analogously to [1]

$$\beta(w, w\_\*) = \frac{\bar{\beta}}{(c + |w - w\_\*|)^{\nu}} \tag{4}$$

where *β*¯ > 0, *ν* > 0 and *c* ≥ 0. According to the above contact rate, agents with similar wealth are more likely to interact. The extrapolation of heterogeneous contact rates have been deeply studied in mathematical epidemiology; see [1,43–47] and the references therein.

Finally, the operators *QH <sup>J</sup>*(*fH*, *fJ*), *H*, *J* ∈ {*<sup>S</sup>*, *I*, *V*, *R*} characterize the evolution of the wealth in each compartment due to wealth exchange activities between agents of the same class, or between agents of different classes *H* and *J*. Their form follows the one originally proposed in the Cordier–Pareschi–Toscani model [4]. An interaction between two individuals in compartment *H* and *J* with wealth pair (*<sup>w</sup>*, *<sup>w</sup>*∗) leads to a wealth pair (*wJH*, *<sup>w</sup>H J*) defined by relations

$$\begin{aligned} w'\_{Hf} &= (1 - \lambda\_H)w + \lambda\_I w\_\* + \eta\_{Hf} w \\ w'\_{]H} &= (1 - \lambda\_I)w\_\* + \lambda\_H w + \eta\_{IH} w\_{\*\*} \end{aligned} \tag{5}$$

with *H*, *J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}. In (5) the constants *λ H*, *λJ* ∈ (0, 1) are exchange parameters defining the saving propensities 1 − *λ H* and 1 − *<sup>λ</sup>J*, i.e., the maximum percentage of money that individuals are willing to employ in a general monetary transaction. Note that the parameters are different in each compartment, underlining the differing behavior of agents in the presence of the pandemic. The choice *λV* > *λS*, for example, reflects the fact that susceptible non-vaccinated agents have reduced action in wealth exchanges due to various governmen<sup>t</sup> restrictions with respect to vaccinated individuals.

Furthermore, *ηJH* ≥ − *λ H*, *ηH J* ≥ − *λJ* are independently centered random variables with the same distribution Θ such that Var(*ηH J*) = Var(*ηH J*) = *<sup>σ</sup>*<sup>2</sup>(*t*). The quantity *σ*<sup>2</sup>(*t*) represents the market risk, which is the same for the whole population and is influenced by the progress of the pandemic. This is in agreemen<sup>t</sup> with market reactions that have been observed during new epidemic waves; see, e.g., ref. [14]. It is convenient to express the operators *QH <sup>J</sup>*(*fH*, *fJ*) in weak form, i.e., the way these operators act on observable quantities [27].

Let *ϕ*(*w*) be a test function and let · denote the expectation with respect to the pair of random variables *ηJH*, *ηH J* in the interaction process (5). Then, for *H*, *J* ∈ {*<sup>S</sup>*, *I*, *V*, *R*} we define the Boltzmann-type bilinear operators as follows

$$\int\_{\mathbb{R}\_+} \varrho(w) Q\_{Hl}(f\_{Hr} f\_l)(w, t) dw = \left\langle \int\_{\mathbb{R}\_+^2} (\varrho(w\_{Hl}') - \varrho(w)) f\_{H}(w, t) f\_l(w\_\*, t) dw \, dw\_\* \right\rangle \tag{6}$$

where (*<sup>w</sup>*, *<sup>w</sup>*∗) → (*wJH*, *<sup>w</sup>H J*) as in (5) and where · denotes the expectation with respect to the independent random variables *ηH J*, *ηH J*.

Binary interactions between individuals (5) reflect the idea that wealth exchanges occur between pairs of agents who invest a fraction of their wealth in the presence of an equivalent good. In each case, such investments involve nondeterministic speculative risks that can provide additional wealth or a loss of wealth. The aggregate behavior of the population is then provided by the operators (6), from which we obtain the emerging macroscopic trends of the binary exchanges considered in each epidemiological compartment.

**Remark 1.** *In the kinetic epidemic model* (2) *the passage from susceptible to vaccinated is governed by a very simple dynamics that does not take into account possible vaccine limitations, as in the first phase of the vaccination campaign. In general, the vaccination rate α may depend on several factors such as the age and work status of individuals and time. It is worthwhile to observe that, in addition to the natural dependency of the recovery rate γI from age [8,22,48], we may also consider wealth-dependent recovery rates to take into account the fact that high wealth can provide access to better hospitals in some health systems, thus ensuring a higher chance of recovery [1]. We point the interested reader to [39] for a more detailed discussion based on the available data.*

#### *2.2. Evolution of Macroscopic Quantities*

In the following, we discuss the evolution of emerging macroscopic quantities from the kinetic model (2). Let *ϕ*(*w*) be a test function. Choosing *ϕ*(*w*) = 1 in (6), we have

$$\sum\_{J \in \{S, I, V, \mathcal{R}\}} \int\_{\mathbb{R}\_+} \varphi(w) Q\_{H \mid \!(f\_{H \prime} f \!\!/)}(w, t) dw = 0,$$

which corresponds to mass conservation, i.e., the conservation of the number of agents. If *ϕ*(*w*) = *w* in (6), we ge<sup>t</sup> the evolution of the average wealth in each compartment, corresponding to the first quantity not conserved in time:

$$\begin{split} \frac{d}{dt} m\_{1,H}(t) &= \frac{1}{H(t)} \sum\_{J \in \{S,I,V,R\}} \int\_{\mathbb{R}^2\_+} \langle w'\_{HI} - w \rangle f\_H(w,t) f\_I(w\_\*,t) dw dw\_\* \\ &= H(t) \sum\_{J \in \{S,I,V,R\}} J(t) (\lambda\_J m\_{1,J}(t) - \lambda\_H m\_{1,H}). \end{split} \tag{7}$$

The total mean wealth is then conserved:

$$\frac{d}{dt}\sum\_{H\in\{S,I,V,R\}}\int\_{\mathbb{R}\_+} w f\_H(w,t)dw = \frac{d}{dt}m\_1 = 0.$$

The evolution of mass fractions can be easily obtained from (2) via direct integration

$$\begin{split} \frac{d}{dt}S(t) &= -\int\_{\mathbb{R}^{2}\_{+}} \beta(w, w\_{\ast}) f\_{\mathbb{S}}(w, t) f\_{\mathbb{I}}(w, t) dw \, dw\_{\ast} - a S(t), \\ \frac{d}{dt} I(t) &= \int\_{\mathbb{R}^{2}\_{+}} \beta(w, w\_{\ast}) f\_{\mathbb{S}}(w, t) f\_{\mathbb{I}}(w, t) dw \, dw\_{\ast} + (1 - \xi) \int\_{\mathbb{R}^{2}\_{+}} \beta(w, w\_{\ast}) f\_{\mathbb{I}}(w, t) f\_{\mathbb{I}}(w, t) dw \, dw\_{\ast} - \gamma\_{\mathbb{I}} I(t), \\ \frac{d}{dt} V(t) &= a S(t) - (1 - \xi) \int\_{\mathbb{R}^{2}\_{+}} \beta(w, w\_{\ast}) f\_{\mathbb{I}}(w, t) f\_{\mathbb{I}}(w, t) dw \, dw\_{\ast}, \\ \frac{d}{dt} R(t) &= \gamma\_{\mathbb{I}} I(t). \end{split} \tag{8}$$

To obtain a closed-form evolution of the macroscopic quantities, we consider a constant rate function, *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) = *β*¯ > 0, obtained from (4) for *ν* = 0, and a constant-in-time market risk *σ*<sup>2</sup>(*t*) = *σ*2. Under these assumptions, thanks tothe mass conservation of Boltzmann-type operators (6), we obtain a classical SIR model with vaccination

$$\begin{aligned} \frac{d}{dt}S(t) &= -\beta S(t)I(t) - aS(t),\\ \frac{d}{dt}I(t) &= \beta S(t)I(t) + (1-\zeta)\tilde{\beta}V(t)I(t) - \gamma\_I I(t),\\ \frac{d}{dt}V(t) &= aS(t) - (1-\zeta)\tilde{\beta}V(t)I(t),\\ \frac{d}{dt}R(t) &= \gamma\_I I(t). \end{aligned} \tag{9}$$

As a consequence, for large times *t* → +<sup>∞</sup>, we have a disease-free equilibrium state, where *I*(*t*) → 0+, *S*(*t*) → 0+, *V*(*t*) → *V*∞ and *R*(*t*) → *R*∞ with *V*∞ + *R*∞ = 1 (see [3]).

The dynamics of mean wealth can be recovered from (7) as follows

$$\begin{aligned} S(t)\frac{d}{dt}m\_{1,S}(t) &= S(t)(\dot{m}\_1(t) - \lambda\_S m\_{1,S}(t)),\\ I(t)\frac{d}{dt}m\_{1,I}(t) &= \vec{\beta}S(t)I(t)(m\_{1,S} - m\_{1,I}) + \vec{\beta}(1 - \vec{\varepsilon})V(t)I(t)(m\_{1,V} - m\_{1,I})\\ &+ I(t)(\dot{m}\_1 - \lambda\_I m\_{1,I}),\\ V(t)\frac{d}{dt}m\_{1,V}(t) &= aS(t)(m\_{1,S} - m\_{1,V}) + V(t)(\dot{m}\_1 - \lambda\_V m\_{1,V}),\\ R(t)\frac{d}{dt}m\_{1,R}(t) &= \gamma\_I I(t)(m\_{1,R}(t) - m\_{1,I}(t)) + R(t)(\dot{m}\_1(t) - \lambda\_R m\_{1,R}(t)),\end{aligned} \tag{10}$$

where we defined the weighted mean wealth as

$$\mathfrak{m}\_1(t) = \sum\_{J \in \{S, I, V, R\}} \lambda\_J m\_{1,J}(t) f(t). \tag{11}$$

Therefore, based on (10), we can observe that the large time behavior of the mean wealth satisfies

$$2\bar{m}\_1^{\infty} - \lambda\_V m\_{1,V}^{\infty} - \lambda\_R m\_{1,R}^{\infty} = 0.$$

Hence, we obtain

$$
\lambda\_V m\_{1,V}^{\infty} = \lambda\_R m\_{1,R'}^{\infty}
$$

together with the constraint *<sup>R</sup>*<sup>∞</sup>*m*<sup>∞</sup>*R*,<sup>1</sup> + *<sup>V</sup>*<sup>∞</sup>*m*<sup>∞</sup>*V*,<sup>1</sup> = *m*, based on the conservation of total mean wealth. Thanks to the latter equalities, we can observe that the asymptotic mean wealth in the compartments of vaccinated and recovered individuals is given by

$$m\_{1,V}^{\infty} = \frac{\lambda\_R}{\lambda\_R V^{\infty} + \lambda\_V R^{\infty}} m\_\prime \qquad m\_{1,R}^{\infty} = \frac{\lambda\_V}{\lambda\_R V^{\infty} + \lambda\_V R^{\infty}} m\_\prime \tag{12}$$

Likewise, we obtain the system for the the second moments

$$\begin{aligned} S(t)\frac{d}{dt}m\_{2,S}(t) &= (\lambda\_S^2 - 2\lambda\_S + \sigma^2)Sm\_{2,S} + S(t)\dot{m}\_2 + 2(1-\lambda\_S)Sm\_{1,S}\dot{m}\_{1,I} \\ I(t)\frac{d}{dt}m\_{2,I}(t) &= \bar{\beta}SI(m\_{2,S} - m\_{2,I}) + (1-\xi)\bar{\beta}VI(m\_{2,V} - m\_{2,I}) \\ &+ (\lambda\_I^2 - 2\lambda\_I + \sigma^2)Im\_{2,I} + Im\_2 + 2(1-\lambda\_I)Im\_{1,I}\dot{m}\_{1,I} \\ V(t)\frac{d}{dt}m\_{2,V}(t) &= aS(m\_{2,S} - m\_{2,V}) + (\lambda\_V^2 - 2\lambda\_V + \sigma^2)Vm\_{2,V} + V\dot{m}\_2 \\ &+ 2(1-\lambda\_V)Vm\_{1,V}\dot{m}\_{1,I} \end{aligned} \tag{13}$$

$$R(t)\frac{d}{dt}m\_{2,R}(t) = (\lambda\_R^2 - 2\lambda\_R + \sigma^2)Rm\_{2,R} + R\dot{m}\_2 + 2(1-\lambda\_R)Rm\_{1,R}\dot{m}\_{1,I}$$

where *m*¯ 1 has been defined in (11) and we have introduced the following notation

$$m\_2(t) = \sum\_{J \in \{S, I, V, R\}} \lambda\_f^2 m\_{2,J}(t) J(t) \dots$$

The evolution of the second moment for the whole system is governed by

$$\frac{d}{dt}\mathfrak{m}\_2(t) = \mathfrak{m}\_2(t) + \sum\_{\substack{I \in \{S,I,V,R\}}} \left( m\_{I,2}(\lambda\_I^2 - 2\lambda\_I + \sigma) + 2(1-\lambda\_I)m\_{I,2}\mathfrak{m}\_1(t) \right) I(t).$$

For large times, the second-order moment for susceptible and infected is such that *<sup>m</sup>*2,*S*, *<sup>m</sup>*2,*I* → 0<sup>+</sup> for *t* → +<sup>∞</sup>. Therefore, *<sup>m</sup>*<sup>∞</sup>2,*V*, *<sup>m</sup>*<sup>∞</sup>2,*R* are solutions to

$$\begin{aligned} (\lambda\_V^2 - 2\lambda\_V + \sigma^2)m\_{2,V}^{\infty} + \vec{m}\_2^{\infty} + (1 - \lambda\_V)m\_{1,V}^{\infty}\vec{m}\_1^{\infty} &= 0, \\ (\lambda\_R^2 - 2\lambda\_R + \sigma^2)m\_{2,R}^{\infty} + \vec{m}\_2^{\infty} + (1 - \lambda\_R)m\_{1,R}^{\infty}\vec{m}\_1^{\infty} &= 0. \end{aligned}$$

from which we ge<sup>t</sup>

$$\begin{split} m\_{2,R}^{\infty} &= \frac{\lambda\_V^2 (1 - \lambda\_V) V^{\infty} m\_{1,V}^{\infty} \bar{m}\_1^{\infty} - A\_V (1 - \lambda\_R) m\_{1,R}^{\infty} \bar{m}\_1^{\infty}}{A\_V (\lambda\_R^2 (1 + R^{\infty}) - 2\lambda\_R + \sigma^2) - \lambda\_V^2 \lambda\_R^2 V^{\infty} R^{\infty}} \\ m\_{2,V}^{\infty} &= \frac{\lambda\_R^2 (1 - \lambda\_R) R^{\infty} m\_{1,R}^{\infty} \bar{m}\_1^{\infty} - A\_R (1 - \lambda\_V) m\_{1,V}^{\infty} \bar{m}\_1^{\infty}}{A\_R (\lambda\_V^2 (1 + V^{\infty}) - 2\lambda\_V + \sigma^2) - \lambda\_R^2 \lambda\_V^2 V^{\infty} R^{\infty}} \end{split}$$

where

$$A\_H = \lambda\_V^2 (1 + H^{\infty}) - 2\lambda\_V + \sigma^2, \qquad H \in \{V, R\}\_{\prime\prime}$$

and *m*¯ ∞1 = *<sup>λ</sup>Vm*<sup>∞</sup>1,*VV*<sup>∞</sup> + *<sup>λ</sup>Rm*<sup>∞</sup>1,*RR*<sup>∞</sup> and *<sup>m</sup>*<sup>∞</sup>1,*V*, *<sup>m</sup>*<sup>∞</sup>1,*R* have been obtained in (12).

**Remark 2.** *In the general case where a non-constant incidence rate β* = *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) *is considered, the macroscopic system of equations is not closed. Depending on the specific choice of β and using the knowledge on the equilibrium states discussed in Section 3.1 it is possible, through the classical hydrodynamic closure of kinetic theory, to derive epidemic models where the dynamics, instead of being homogeneous as in classical compartmental modeling, is influenced by the heterogeneous wealth status of individuals. We refer to [8,38] for examples in this direction.*

#### **3. Properties of the Kinetic Model**

In this section we study the mathematical model (2) from an analytical point of view, by proving the well-posedness and convergence to equilibrium of the solution. To this end, we made suitable simplification assumptions on the contact rate by restricting to the case *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) = *β*¯. We resort to classical mathematical approaches for kinetic equations to characterize the trend to equilibrium [1,27]. In particular, taking into account methods

for nonconservative systems—see, e.g., ref. [49]—we provide an existence and uniqueness result. Given a function *f*(*w*) ∈ *<sup>L</sup>*<sup>1</sup>(R+), we define its Fourier transform as follows

$$\widehat{f}(z) = \int\_{\mathbb{R}} e^{-iwz} f(w) dw.$$

According to the above assumption regarding the contact rate, we rewrite (2) in weak form:

*∂t* R+ *ϕ*(*w*)*fS*(*<sup>w</sup>*, *t*)*dw* = −*β*¯*I*(*t*) R+ *ϕ*(*w*)*fS*(*<sup>w</sup>*, *t*)*dw* − *α* R+ *ϕ*(*w*)*fS*(*<sup>w</sup>*, *t*)*dw* + ∑ *J*∈{*<sup>S</sup>*,*I*,*V*,*<sup>R</sup>*} R+ *ϕ*(*w*)*QS <sup>J</sup>*(*fS*, *fJ*)(*<sup>w</sup>*, *<sup>t</sup>*)*dw*, *∂t* R+ *ϕ*(*w*)*fI*(*<sup>w</sup>*, *t*)*dw* = *β*¯*I*(*t*) R+ *ϕ*(*w*)*fS*(*<sup>w</sup>*, *t*)*dw* + (1 − *ζ*)*β*¯*I*(*t*) R+ *ϕ*(*w*)*fV*(*<sup>w</sup>*, *t*)*dw* − *γI* R+ *ϕ*(*w*)*fI*(*<sup>w</sup>*, *t*)*dw* + ∑ *J*∈{*<sup>S</sup>*,*I*,*V*,*<sup>R</sup>*} R+ *ϕ*(*w*)*QI <sup>J</sup>*(*fI*, *fJ*)(*<sup>w</sup>*, *<sup>t</sup>*)*dw*, *∂t* R+ *ϕ*(*w*)*fV*(*<sup>w</sup>*, *t*)*dw* = *α* R+ *ϕ*(*w*)*fS*(*<sup>w</sup>*, *t*)*dw* − (1 − *ζ*)*β*¯*I*(*t*) R+ *ϕ*(*w*)*fV*(*<sup>w</sup>*, *t*)*dw* + ∑ *J*∈{*<sup>S</sup>*,*I*,*V*,*<sup>R</sup>*} R+ *ϕ*(*w*)*QV <sup>J</sup>*(*fV*, *fJ*)(*<sup>w</sup>*, *<sup>t</sup>*)*dw*, *∂t* R+ *ϕ*(*w*)*fR*(*<sup>w</sup>*, *t*)*dw* = *γI* R+ *ϕ*(*w*)*fI*(*<sup>w</sup>*, *t*)*dw* + ∑ *J*∈{*<sup>S</sup>*,*I*,*V*,*<sup>R</sup>*} R+ *ϕ*(*w*)*QR <sup>J</sup>*(*fR*, *fJ*)(*<sup>w</sup>*, *<sup>t</sup>*)*dw*. (14)

Hence, we consider *ϕ*(*w*) = *e*<sup>−</sup>*izw* in (14) to ge<sup>t</sup>

$$\begin{split} \partial\_{t} f\_{S}(z,t) &= -\not{\partial} I(t) f\_{S}(z,t) - a f\_{S}(z,t) + \sum\_{I \in \{S,I,V,R\}} \dot{Q}\_{Sf}(\hat{f}\_{S},f\_{I})(z,t), \\ \partial\_{t} \hat{f}\_{I}(z,t) &= \not{\partial} I(t) \hat{f}\_{S}(z,t) + (1-\zeta) \not{\partial}\_{I}^{2} f\_{I}(z,t) \hat{f}\_{V}(z,t) - \gamma\_{I} \not{f}\_{I}(z,t) + \sum\_{I \in \{S,I,V,R\}} \dot{Q}\_{II}(\hat{f}\_{I},\hat{f}\_{I})(z,t), \\ \partial\_{t} \hat{f}\_{V}(z,t) &= a \hat{f}\_{S}(z,t) - (1-\zeta) \not{\partial}\_{I}^{2} f\_{I}(z,t) \hat{f}\_{V}(z,t) + \sum\_{I \in \{S,I,V,R\}} \dot{Q}\_{VI}(\hat{f}\_{V},\hat{f}\_{I})(z,t), \\ \partial\_{t} \hat{f}\_{R}(z,t) &= \gamma\_{I} \dot{f}\_{I}(z,t) + \sum\_{I \in \{S,I,V,R\}} \dot{Q}\_{Rf}(\hat{f}\_{R},\hat{f}\_{I})(z,t). \end{split} \tag{15}$$

Similarly to [1] the operators *Q* ˆ *H J*( ˆ*fH*, ˆ*fJ*)(*<sup>z</sup>*, *t*) may be rewritten as follows

$$\int\_{\mathbb{R}\_+} e^{-iwz} Q\_{H\{}}(f\_{H\nu}f\_{\}) dw = \langle \hat{f}\_H(A\_{H\{}^c}z, t) \rangle \hat{f}\_{\{}^c}(\lambda\_{\} z, t) - f(t) \hat{f}\_H(z, t) \rangle$$

where

$$A\_{HI} = 1 - \lambda\_H + \eta\_{HI}.$$

We assume that the parameters of the trading activity satisfy the condition

$$\nu = \max\_{H,l \in \{S,I,V,\mathcal{R}\}} \left[ \lambda\_f^2 + \langle A\_{HI}^2 \rangle \right] < 1. \tag{16}$$

Let P*s*(R+) be the set of probability measures *f*(*w*) with bounded *<sup>s</sup>*−moment, and, for any pair of densities *f* and *g* in P*s*(R+), let us consider the class of metrics *ds* defined by

*ds*(*f* , *g*) = sup *z*∈R | ˆ*f*(*z*) − *g*<sup>ˆ</sup>(*z*)| |*z*|*<sup>s</sup>* , (17)

where ˆ *f* and *g*ˆ denote the Fourier transforms of *f* and *g*. Then, the distance (17) is welldefined and finite for any pair of probability measures with equal moments up to order

[*s*] (where [*s*] denotes the integer part of *s*), if *s* is a real number or up to *s* − 1, if *s* is an integer [27].

Inequality (16), combined with a Fourier-based distance, allows one to obtain an exponential convergence to equilibrium for system (2). This condition is verified whenever

$$
\sigma^2 < 2 \min\_{J \in \{S, I, V, R\}} \lambda\_J (1 - \lambda\_J)\_J
$$

namely, when the market risk is not too grea<sup>t</sup> in relation to the saving propensities. To study the large-time behavior of the solution to systems such as (15) we follow [1,27]. Then, we have the following result

**Theorem 1.** *Let fJ*(*<sup>w</sup>*, *t*) *and gJ*(*<sup>w</sup>*, *t*)*, J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}*, be two solutions of the kynetic system* (2)*, corresponding to the initial values fJ*(*<sup>w</sup>*, 0) *and gJ*(*<sup>w</sup>*, 0) *such that <sup>d</sup>*2(*fJ*(*<sup>w</sup>*, <sup>0</sup>), *gJ*(*<sup>w</sup>*, <sup>0</sup>))*, J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}*, is finite. Then, if condition* (16) *holds, the Fourier-based distance <sup>d</sup>*2(*fJ*(*<sup>w</sup>*, *t*), *gJ*(*<sup>w</sup>*, *t*)) *decays exponentially in time toward zero and the following holds:*

$$\sum\_{\{f\in\{S,I,V,R\}}} d\_2(f\_I(w,t),g\_I(w,t)) < \sum\_{f\in\{S,I,V,R\}} d\_2(f\_I(w,0),g\_I(w,0)) \exp\{- (1-\nu)t\}. \tag{18}$$

The previous result and the Equation (18) give us the contractivity of the system in the *d*2 metric, which will be the essential to prove the existence theorem. Theorem 1 allows us to further investigate the properties of the steady state *f* ∞*J* (*w*), *J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}.

In order to obtain an existence result we need to introduce a subset of P2(R)

$$\mathcal{D}\_{m\_1, m\_2} := \left\{ F \in \mathcal{P}\_2(\mathbb{R}) : \int\_{\mathbb{R}} v dF(v) = m\_1, \int\_{\mathbb{R}} v^2 dF(v) = m\_2 \right\}. \tag{19}$$

Following [49], it is possible to prove that <sup>D</sup>*<sup>m</sup>*1,*m*2 is a metric Banach space with the *d*2(·, ·) metric. Now, we define

$$\mathcal{D}^{\infty} := \mathcal{D}\_{\mathfrak{m}\_{V,1}^{\infty}, \mathfrak{m}\_{V,2}^{\infty}} \times \mathcal{D}\_{\mathfrak{m}\_{R,1}^{\infty}, \mathfrak{m}\_{R,2}^{\infty}}$$

as the product space of two sets such as (19), where the momenta are those of the steady states for the relative distributions *fJ*(*w*), for *J* ∈ {*<sup>V</sup>*, *R*} (we are only considering these two classes since for large time *I*, *S* → 0<sup>+</sup>). We also recall a variant of the metric used in Theorem 1

$$\overline{d}\_2(f, \emptyset) := \sum\_{J \in \{V, R\}} d\_2(f\_I(w, t), \emptyset\_J(w, t)). \tag{20}$$

Now, we are able to prove the following theorem.

**Theorem 2.** *If the initial value f*0(*w*) = *f*(*<sup>w</sup>*, 0) ∈ D∞ *and condition* (16) *holds, then the system*

$$\begin{aligned} \partial\_t f\_V(w, t) &= \sum\_{J \in \{V, R\}} Q\_{Vf}(f\_{V, \prime} f\_I)(w, t), \\ \partial\_t f\_R(w, t) &= \sum\_{J \in \{V, R\}} Q\_{Rf}(f\_{R\prime} f\_I)(w, t), \end{aligned} \tag{21}$$

*has a unique steady state f* <sup>∞</sup>(*w*)*, and it also belongs to* D<sup>∞</sup>*.*

**Proof.** Let us consider the flow map

$$T\_{\rm l}: \left(\mathcal{D}^{\infty}, \bar{d}\_2\right) \to \left(\mathcal{D}^{\infty}, \bar{d}\_2\right) \tag{22}$$

which , for any time *t* > 0, is given by *Tt*(*f*0(*w*)) = *f*(*t*)=(*fV*(*<sup>w</sup>*, *t*), *fR*(*<sup>w</sup>*, *<sup>t</sup>*)), where *f*(*t*) is the solution of (21) at time *t* with *f*(*<sup>w</sup>*, 0) = *f*0(*w*) ∈ D<sup>∞</sup>. Thanks to (18) we have

$$d\_2\left(T\_t(f\_0(w)), \ T\_t(\mathcal{g}\_0(w))\right) < d\_2\left(f\_0(w), \ \mathcal{g}\_0(w)\right) \exp\{- (1 - \nu)t\}$$

which is a strict contraction for (22) with constant exp{−(<sup>1</sup> − *ν*)*t*} < 1. Now, it is easy to see that -D<sup>∞</sup>, *<sup>d</sup>*2 is a Banach space and therefore the Banach fixed-point theorem ensures the existence and uniqueness for the steady state in D<sup>∞</sup>.

**Remark 3.** *Similar results may be obtained in the more realistic case β*(*<sup>w</sup>*, *<sup>w</sup>*∗) = *β*(*w* − *<sup>w</sup>*∗) *since the transmission operator <sup>K</sup>*(·, ·) *defined in* (3) *possesses, in this case, a convolution structure, which naturally converts into a product in the Fourier space. We omit the details.*

#### *3.1. Fokker–Planck Scaling and Steady States*

In the general case, it is difficult to compute analytically the large-time behaviour of the compartmental kinetic system (2). A deeper insight into the steady states can be obtained through the so-called quasi-invariant limit procedure [1,4,27]. The goal is to derive a simplified Fokker–Planck model in which the study of the asymptotic properties is much easier. It is worth mentioning that this approach is inspired by the so-called grazing collision limit of the Boltzmann equation; see [50,51].

The driving idea is to scale interactions and trading frequency at the same time. As a consequence, the equilibrium of the wealth distribution is reached more quickly with respect to the time scale of the epidemic. Hence, given 1 we introduce the following scaling

$$\begin{aligned} \lambda\_S &\rightarrow \epsilon \lambda\_{S\prime} & \lambda\_I &\rightarrow \epsilon \lambda\_{I\prime} & \lambda\_V &\rightarrow \epsilon \lambda\_{V\prime} & \lambda\_R &\rightarrow \epsilon \lambda\_{R\prime} \\ \sigma^2 &\rightarrow \epsilon \sigma^2, & \beta(w, w\_\*) &\rightarrow \epsilon \beta(w, w\_\*), & \gamma\_I &\rightarrow \epsilon \gamma\_{I\prime} \end{aligned} \tag{23}$$

together with the time scaling *t* → *t*/. We denote as *QH <sup>J</sup>*(·, ·), *H*, *J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}, the scaled interaction terms. Using a Taylor expansion for small values of , we ge<sup>t</sup> [1]

$$\begin{aligned} &\frac{1}{2}\int\_{\mathbb{R}\_+} Q\_{H\_I}^{\varepsilon}(f\_{H\_I}f\_I)(w,t)\,\eta(w)dw \\ &= \int\_{\mathbb{R}\_+} \left\{-q^{\prime}(w)(w\lambda\_H f - m\_{1,l}\lambda\_I) + \frac{\sigma^2}{2}q^{\prime\prime}(w)w^2I(t)\right\}f\_H(w,t)dw + O(\varepsilon). \end{aligned}$$

Integrating back by parts, in the limit → 0, we obtain the system of Fokker–Planck equations

$$\begin{aligned} \frac{\partial f\_S(w,t)}{\partial t} &= -K(f\_S, f\_I)(w,t) - a f\_S(w,t) + \frac{\partial}{\partial w} \{ [w\lambda\_S - \dot{m}(t)] f\_S(w,t) \} \\ &+ \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} (w^2 f\_S(w,t)),\\ \frac{\partial f\_I(w,t)}{\partial t} &= K(f\_S, f\_I)(w,t) + (1-\zeta)K(f\_V, f\_I)(w,t) - \gamma\_I f\_I(w,t) \\ &+ \frac{\partial}{\partial w} \{ [w\lambda\_I - \dot{m}(t)] f\_I(w,t) \} + \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} (w^2 f\_I(w,t)),\\ \frac{\partial f\_V(w,t)}{\partial t} &= a f\_S(w,t) - (1-\zeta)K(f\_V, f\_I)(w,t) + \frac{\partial}{\partial w} \{ [w\lambda\_V - \dot{m}(t)] f\_V(w,t) \} \\ &+ \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} (w^2 f\_V(w,t)),\\ \frac{\partial f\_R(w,t)}{\partial t} &= \gamma\_I(w,t) + \frac{\partial}{\partial w} \{ [w\lambda\_R - \dot{m}(t)] f\_R(w,t) \} + \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} (w^2 f\_R(w,t)), \end{aligned} \tag{24}$$

where *m*¯ has been defined in (11). The above Fokker–Planck system is complemented with the following boundary conditions

$$\frac{\partial}{\partial w} [w^2 g\_J(w, t)]|\_{w=0} = 0 \qquad [w\lambda\_J - \overline{m}] g\_J + \frac{\sigma}{2} \frac{\partial}{\partial w} (w^2 g\_J) \Big|\_{w=0} = 0.$$

We can verify under suitable assumptions that the Fokker–Planck system (24) possesses an explicitly computable steady state [52]. Let us consider the case of a constant contact rate, i.e., *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) = *β*¯. Since for large times *S*, *I* → 0<sup>+</sup> we find that the stationary states *f* ∞*V* (*w*) and *f* ∞*R*(*w*) solve the following equations:

$$\begin{aligned} \lambda\_V \frac{\partial}{\partial w} \left[ (w - m\_V^{\infty}) f\_V^{\infty} (w) \right] + \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} [w^2 f\_V^{\infty} (w)] &= 0, \\ \lambda\_R \frac{\partial}{\partial w} \left[ (w - m\_R^{\infty}) f\_R^{\infty} (w) \right] + \frac{\sigma^2}{2} \frac{\partial^2}{\partial w^2} [w^2 f\_R^{\infty} (w)] &= 0. \end{aligned}$$

Based on the above equalities, we find that the two steady states are inverse Gamma densities

$$f\_V^{\infty}(w) = V^{\infty} \frac{\kappa^{\mu\_V}}{\Gamma(\mu\_V)} \frac{e^{-\frac{\mathcal{E}}{w}}}{w^{1+\mu\_V}} \qquad \qquad f\_R^{\infty}(w) = R^{\infty} \frac{\kappa^{\mu\_R}}{\Gamma(\mu\_R)} \frac{e^{-\frac{\mathcal{E}}{w}}}{w^{1+\mu\_R}} \tag{25}$$

with Pareto indices defined as follows

$$\mu\_V = 1 + 2\frac{\lambda\_V}{\sigma^2}, \qquad \mu\_R = 1 + 2\frac{\lambda\_R}{\sigma^2}.$$

$$\kappa = (\mu\_V - 1)m\_V^{\infty} = (\mu\_R - 1)m\_R^{\infty} = \frac{2\lambda\_R\lambda\_V}{\sigma^2(\lambda\_R V^{\infty} + \lambda\_V R^{\infty})}m.$$

Consequently, the global steady state is a mixture of the inverse Gamma distribution

$$f^{\infty}(w) = f\_V^{\infty}(w) + f\_R^{\infty}(w),\tag{26}$$

which may present a bimodal shape with a different intensity.The formation of two peaks at the equilibrium is due to the fact that we have two different maxima corresponding to the points

$$
\overline{w}\_V = \frac{\kappa}{\mu\_V + 1} = \frac{\lambda\_R \lambda\_V}{(\lambda\_V + \sigma)(\lambda\_R V^{\infty} + \lambda\_V R^{\infty})} m\_\prime \tag{27}
$$

$$\overline{w}\_{R} = \frac{\kappa}{\mu\_{R} + 1} = \frac{\lambda\_{R}\lambda\_{V}}{(\lambda\_{R} + \sigma)(\lambda\_{R}V^{\infty} + \lambda\_{V}R^{\infty})} \, m\_{\prime} \tag{28}$$

for the vaccinated and for the recovered wealth distributions, respectively. In the next section we report on the resulting profiles for different choices of *λV*, *λR*, *σ* and *V*<sup>∞</sup>, *R*<sup>∞</sup>.

**Remark 4.** *The emergence of a multimodal equilibrium wealth distribution has been classically linked to the appearance of new inequalities in highly stressed societies; see, e.g., [15,35,53]. In these cases, the economic segregation of part of the society leads to the pauperization of substantial layers of the middle class. In the present case, the different economic impact played by agents in each compartment is capable of shaping the wealth distribution towards a bimodal distribution. Indeed, the trading propensities modeling personal responses to the economic scenario can be substantially modified by the progression of the epidemic and the vaccine efficacy.*

#### **4. Numerical Results**

In this section we study the impact of vaccination on the equilibrium of the kinetic system through several numerical simulations. This allows us to show the model's ability to describe different situations of wealth distribution in the presence of epidemic dynamics. In particular, we will adopt standard direct simulation Monte Carlo methods

to simulate the system of kinetic Equation (2); see [27] and the references therein. In all the subsequent tests we will consider *N* = 10<sup>5</sup> agents and the densities are reconstructed through standard histograms.

In the first test, we verify numerically the convergence of the solution to the kinetic system (2) to the solution of the Fokker–Planck system (24) under the scaling (23). Then, we study the emergence of wealth inequalities, measured through the Gini index, in relation to the effectiveness of the vaccine. These results are obtained both in the case of a constant market risk variance *σ*<sup>2</sup> and in the case of a variance that depends on the current epidemic situation. Lastly, we introduce the possibility that the effectiveness of the vaccine is also affected by the number of positive cases. This situation mimics the realistic case of the diffusion of viral variants for which an up-to-date vaccine may be not immediately available.

#### *4.1. Test 1: Long-Time Behavior and Convergence to Equilibrium*

In this test, we want to observe the convergence of the numerical solution of the kinetic system (2) to the one of the Fokker–Planck system (24) in the quasi-invariant limit introduced in Section 3.1. We consider the simplified case where *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) = *β*¯ = 0.2, *γI* = 1/12 and *ζ* = 0.9, for which we obtained the steady distributions in (25). These values are representative of realistic dynamics during the beginning of the COVID-19 pandemic; see, e.g., [6–8,39,54].

At time *t* = 0 we consider an inverse Gamma distribution

$$f(w,0) = \frac{(\mu - 1)^{\mu}}{\Gamma(\mu)} \frac{\exp\left(-\frac{\mu - 1}{\mu}\right)}{w^{1 + \mu}},\tag{29}$$

where <sup>Γ</sup>(·) is the Gamma function and *μ* = 10. The distributions of the epidemic compartments are

$$f\_S(w,0) = \rho\_S f(w), \quad f\_l(w,0) = \rho\_l f(w), \quad f\_V(w,0) = \rho\_V f(w), \quad f\_R(w,0) = \rho\_R f(w), \tag{30}$$

where the mass fractions are *ρI* = 7.5 × <sup>10</sup>−3, *ρV* = 0, *ρR* = 4 × 10−<sup>2</sup> and *ρS* = 1 − (*ρI* + *ρV* + *ρR*). Furthermore, we consider the value *σ*<sup>2</sup> = 0.02 for the market risk. In Figure 1 we show the numerical solution at time *T* = 300 of (2) in the scaling regime (23) with = 1, 0.5, 10−3.

In particular, provided an epidemic dynamics such that *V* ∞ = 0.51 and *R* ∞ = 0.49, we give numerical evidence of the aforementioned convergence in two regimes expressing increasing safeguard thresholds 1 − *<sup>λ</sup>J*, *J* ∈ {*<sup>S</sup>*, *I*, *V*, *<sup>R</sup>*}, for non-vaccinated agents

 $\lambda(i) \quad \lambda\_S = 0.15$ ,  $\lambda\_I = 0.10$ ,  $\lambda\_V = 0.30$ ,  $\lambda\_R = 0.20$ 

(*ii*) *λS* = 0.10, *λI* = 0.05, *λV* = 0.30, *λR* = 0.15

where the same values of *V* ∞ and *R* ∞ are unchanged. In particular, we assume that recovered individuals are characterized by a greater safeguard parameter. This is coherent with the possibility of reinfection, which will be investigated in the last numerical test.

We observe that, if 1, the Fokker–Planck asymptotic distribution is a consistent approximation of the equilibrium distribution of the Boltzmann-type model. In both cases, the global distribution is a mixture of inverse Gamma densities and in the righ-hand plot depicted in Figure 1, we can clearly observe a bimodal shape for the wealth distribution. To highlight this, we have drawn the maximum points of the distributions *f* ∞ *V* , *f* ∞ *R* , which are at *wV*, *wR*, defined in (27) and (28).

**Figure 1.** Test 1. Comparison of the wealth distributions at the end of the epidemic for the kinetic system (2) with the explicit Fokker–Planck asymptotics (26) with scaling parameters = 1, 12 , 10−3. (**Left**) *λS* = 0.15, *λI* = 0.10, *λV* = 0.30, *λR* = 0.20. (**Right**) *λS* = 0.10, *λI* = 0.05, *λV* = 0.30 *λR* = 0.15. In both cases we fixed *β* ¯ = 0.2, *γI* = 1/12, *α* = 0.005, *ζ* = 0.9 and *σ*<sup>2</sup> = 0.02.

#### *4.2. Test 2: Wealth Inequalities and Vaccination Campaign*

In the second test case we analyze the emergence of wealth inequalities through the computation of the Gini index. In particular, we concentrated on the effects linked to the outbreak of the infection and on the impact of an effective vaccination campaign.

We fixed the epidemic parameters as follows: *β* ¯ = 0.15, *γI* = 1/12 and a vaccination rate of *α* = 10−2. Furthermore, we considered two different vaccine efficacies *ζ* = 0.95, corresponding to a high efficacy of the vaccine, and *ζ* = 0.55 corresponding to a low efficacy of the vaccine. Since we are interested in the behavior of the system up to the conclusion of the epidemic phenomenon, the final time was fixed as *T* = 810, corresponding to a wide time-span. We kept the same values for the saving propensities and market risk as those defined for Section 4.1. Hence, we considered initial wealth distributions as in (29) and mass fractions as in (30), with *ρI* = 7 × <sup>10</sup>−3, *ρV* = 0, *ρR* = 4 × 10−<sup>2</sup> and *ρS* = 1 − (*ρI* + *ρV* + *ρR*). The scaling coefficient was = 5 × 10−2. The resulting epidemic dynamic is reported in Figure 2.

**Figure 2.** Test 2. Evolution of the epidemic dynamics from (9) for the choice of parameters *β* ¯ = 0.15, *γI* = 1/12, *α* = 0.01 and *ζ* = 0.95 (**left**), *ζ* = 0.55 (**right**).

We evaluated the Gini coefficient of the emerging equilibrium distributions. The Gini index is commonly computed from the Lorenz curve

$$L(F(w)) = \int\_0^w f^\infty(w\_\*) w\_\* dw\_\*,$$

where *<sup>F</sup>*(*w*) = *w*0 *f* <sup>∞</sup>(*<sup>w</sup>*∗)*dw*∗ and is defined as follows

$$G\_1 = 1 - 2\int\_0^1 L(\mathbf{x})d\mathbf{x}.$$

This index should be understood as a measure of a country's wealth discrepancy and it varies in [0, 1], where in the case *G*1 = 0 the country is in a situation of perfect equality, whereas *G*1 = 1 indicates complete inequality. A reasonable value for this parameters is in the range [0.2, 0.5] for most Western economies [36].

In Figure 3 we show the evolution of the Gini index with the parameters described above. We may observe that the epidemic peak leads to an increasing of inequalities that is then absorbed for later times in relation to the efficacy of the vaccine. Consequently, only when the vaccine is made available to the majority of the population does it actually contribute to reducing inequalities; otherwise, it may have the opposite effect. This reminds us of how, on a global level, the importance of making vaccines available to all countries should be seen not only in terms of epidemics, but also in terms of reducing economic inequalities. In all the considered cases, in the long term, the Gini index decreases thanks to the vaccine.

**Figure 3.** Test 2. Evolution of Gini index under the epidemic dynamics described in Figure 2 and for the choice of parameters *λS* = 0.10, *λI* = 0.07, *λV* = 0.30, *λR* = 0.15. Two vaccine efficacies were considered: 95% (green) and 55% (red). In both cases we considered *σ*<sup>2</sup> = 0.02.

Next, we consider the case where the market risk is related to the behavior of the epidemic's spread and where there is a linear relation between the market risk and the number of people infected. The introduction of a time-dependent market risk *σ*<sup>2</sup>(*t*) mimics an instantaneous influence of the pandemic on the volatility of a market economy, as is often observed. Therefore, we consider the following:

$$
\sigma^2(t) = \sigma\_0^2 (1 + \mu I(t)) \tag{31}
$$

where *μ* > 0 expresses the effective influence of the epidemic dynamics on the market volatility and *σ*20> 0 is an ineradicable baseline risk.

 In the following, we choose *μ* = 50 and *σ*20 = 0.02. In Figure 4 we represent the evolution of *σ*<sup>2</sup>(*t*) in the presence of an epidemic characterized by *β*¯ = 0.15, *γI* = 1/10. Furthermore, we compare the Gini index in the presence of two effectiveness rates of the vaccine, i.e., *ζ* = 0.95 and *ζ* = 0.55. We may easily observe how an increasing variability leads to a worsening of the Gini index and, therefore, of the inequalities. The long-term behavior of the Gini index depends, as before, on the vaccine efficacy *ζ* such that low efficacy leads to increasing inequalities in the long term. This is due to the fact that as *t* → +∞ we have *I* → 0<sup>+</sup> and then *σ*<sup>2</sup>(*t*) → *σ*20 .

**Figure 4.** Test 2. (**Left**) evolution of the market risk *σ*<sup>2</sup>(*t*) as defined in (31) with *μ* = 50 and *σ*20 = 0.02 in case of two different vaccine efficacies. (**Right**) evolution of Gini index under the epidemic dynamics described in Figure 2 and epidemic-dependent market risk parameter (31).

Finally, in Figure 5 we present the evolution of the full kinetic density solution to (2) in the scaling = 5 × 10−<sup>2</sup> in the presence of fixed market risk *σ*<sup>2</sup> or with the epidemicdependent *σ*<sup>2</sup>(*t*) discussed in (31).

#### *4.3. Nonlinear Incidence Rate and Time-Varying Vaccine Efficacy*

In this last test case, to model different frequencies of interactions between agents that belong to the same wealth class, we introduce a wealth-dependent contact rate *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) of the form ¯

$$\beta(w, w\_\*) = \frac{\beta}{(c + |w - w\_\*|)^{\nu}},\tag{32}$$

where *β* ¯ , *c*, *ν* > 0. We have depicted the above contact rate in Figure 6.

We also introduce a time-dependent efficacy of the vaccine *ζ* of the form

$$\mathcal{L}(t) = \mathcal{\zeta}\_0 - \Psi \int\_0^t \int\_{\mathbb{R}\_+} f\_I(w, t) dw ds = \mathcal{\zeta}\_0 - \Psi \int\_0^t I(s) ds,\tag{33}$$

with *ζ*0 ∈ [0, 1] indicating the initial efficacy of the vaccine and 0 < *ψ* ≤ *ζ*0. This timedependence in vaccine coverage describes, in a simplified way, the fact that with more infected individuals it is more likely to encounter mutations of the original virus, for which the vaccine is less effective. In the following, we compare the evolution of the wealth inequalities in the presence of two different values *ζ*0.

**Figure 5.** Test 2. Time evolution of the wealth distribution of the kinetic model (2) in the scaling = 5 × 10−<sup>2</sup> with vaccine efficacy *ζ* = 0.55 (**left** column) or *ζ* = 0.95 (**right** column) and with constant market risk *σ*<sup>2</sup> = 0.02 (top row) or *<sup>σ</sup>*<sup>2</sup>(*t*), defined in (31) with *μ* = 50. In all the evolutions we considered *λS* = 0.10, *λI* = 0.07, *λV* = 0.30 and *λR* = 0.15. The initial distribution was defined in (29) and (30). In the left image, we can observe the evolution of the wealth distribution for the kinetic model (2) in the scaling parameter = 5 × 10−<sup>2</sup> with *ζ* = 0.95, whereas, in the right image we have the comparison between the behaviors of the Gini index with vaccine effectiveness, equal to 95% (green line) and 65% (red line). In both images we considered a variable market risk (31) with *σ*20= 0.02 and *μ* = 50 and *λS* = 0.10, *λI*= 0.07, *λV* = 0.30 and *λR* = 0.15.

Furthermore, to make the modeling more realistic, we assume the loss of immunity of the agents in the compartment *R*. To this end, we have to modify the first and last equations of the model (2) as follows

$$\begin{aligned} \partial\_t f\_{\mathcal{S}}(w,t) &= -K(f\_{\mathcal{S}},f\_{\mathcal{I}})(w,t) - af\_{\mathcal{S}}(w,t) + \gamma\_{\mathcal{R}}f\_{\mathcal{R}}(w,t) + \sum\_{J \in \{S,I,V,\mathcal{R}\}} Q\_{\mathcal{S}f}(f\_{\mathcal{S}},f\_{\mathcal{I}})(w,t) \\ \partial\_t f\_{\mathcal{R}}(w,t) &= \gamma\_{\mathcal{I}}f\_{\mathcal{I}}(w,t) - \gamma\_{\mathcal{R}}f\_{\mathcal{R}}(w,t) + \sum\_{J \in \{S,I,V,\mathcal{R}\}} Q\_{\mathcal{R}f}(f\_{\mathcal{R}},f\_{\mathcal{I}})(w,t), \end{aligned} \tag{34}$$

where *γR* ≥ 0 is the rate expressing the loss of immunity of recovered agents. Note that this latter assumption substantially changes the epidemic dynamics, since asymptotically, instead of a disease-free scenario, we have the emergence of endemic states [3].

**Figure 6.** Test 3. Wealth-dependent contact rate *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) of the form (32) with *β*¯ = 8, *c* = 7, *ν* = 2.

#### 4.3.1. Test 3A: *γR* = 0

First, we consider model (2) without the modified relations (34) (or equivalently, in the absence of reinfection, i.e., *γR* = 0) and, as before, a fixed recovery rate *γI* = 1/12 and a vaccination rate *α* = 0.005 with the same initial masses as those defined in (30). Furthermore, we fixed *ψ* = 0.005. In Figure 7, in the top row, we show the evolution for the fractions of the population in the case of *ζ*0 = 0.95 (left) and *ζ*0 = 0.55 (right). We may observe how a variable efficacy of the vaccine, affected by epidemic peaks, may strongly shape the immunity of the population, even in the presence of an initial high efficacy. Interestingly, in this latter case, a variable efficacy leads to the emergence of secondary peaks of infection. This is due to the presence of a smaller number of recovered persons who, unlike vaccinated people, maintain immunity.

In Figure 7, in the bottom-left row, we can observe the evolution of the resulting vaccine efficacy for *ζ*0 = 0.95, *ζ*0 = 0.55 and *ψ* = 0.005. The vaccine efficacy is degraded by the epidemic dynamics due to the increasing of the infected compartment, with a slower efficacy decay for high initial *ζ*0.

For the same choice of coefficient, in the bottom-right plot of Figure 7, we show the evolution of the Gini coefficient in the case of variable efficacy as (33). With respect to a vaccine with constant efficacy, the efficacy decay forces the emergence of sharper inequalities, which is well evidenced by the evolution of the Gini coefficient.

**Figure 7.** *Cont*.

**Figure 7.** Test 3A. Top row: epidemic dynamics with wealth-dependent *β*(*<sup>w</sup>*, *<sup>w</sup>*∗), defined in (32) with *β* ¯ = 8, *c* = 7, *ν* = 2, *γI* = 1/12, *α* = 0.005 and variable *ζ* as in (33) with *ψ* = 0.005. We considered *ζ*0 = 0.95 (**left**) and *ζ*0 = 0.55 (**right**). The initial distribution is (29) with mass fractions (30). Bottom row: decline in vaccine efficacy due to the presence of a high number of infective people (**left**) and the evolution of the Gini index (**right**) for a variable infection rate *β*(*<sup>w</sup>*, *<sup>w</sup>*∗) as in (32) and vaccine effectiveness *ζ*(*t*) as in (33). We considered *λS* = 0.10, *λI* = 0.07, *λV* = 0.25, *λR* = 0.15 and *β*¯ = 8, *c* = 7, *ν* = 2 and *ψ* = 0.005.

#### 4.3.2. Test 3B: *γR* > 0

Finally, we consider model (2) including the modified Equation (34), with a reinfection period of 180 days, i.e., *γR* = 1/180 and, as before, a fixed recovery rate *γI* = 1/12 and vaccination rate *α* = 0.005 with the same initial masses as those defined in (30). In the first row of Figure 8, we present two epidemic dynamics with nonlinear contact rates (32) and the time-dependent efficacy *ζ*(*t*) defined in (33) with *ψ* = 1.5 × 10−4. In the left plot, we present the case of strong initial vaccine efficacy *ζ*0 = 0.95 and in the right plot the case of mild initial vaccine efficacy *ζ*0 = 0.45. The macroscopic dynamics present an endemic equilibrium due to the presence of the reinfection rate *γR*. Furthermore, in contrast to the previous case, in the case of reduced initial efficacy of the vaccine, a second infection wave is seen to emerge.

**Figure 8.** *Cont*.

**Figure 8.** Test 3B. Top row: epidemic dynamics with wealth-dependent *β*(*<sup>w</sup>*, *<sup>w</sup>*∗), defined in (32) with *β* ¯ = 8, *c* = 7, *ν* = 2, *γI* = 1/12, *γR* = 1/180, *α* = 0.005 and variable *ζ* as in (33) with *ψ* = 1.5 × 10−4. We considered *ζ*0 = 0.95 (**left**) and *ζ*0 = 0.55 (**right**). The initial distribution is (29) with mass fractions (30). Bottom row: decline in vaccine efficacy due to the presence of a high number of infected people (**left**) and evolution of the Gini index (**right**). We considered *λS* = 0.10, *λI* = 0.07, *λV* = 0.25, *λR* = 0.15 and *β* ¯ = 8, *c* = 7 and *ν* = 2.

Looking at the bottom-left plot, we can observe that, in the present regime of parameters, a strong initial vaccine efficacy is robust with respect to the efficacy decay due to epidemic waves. On the other hand, mild initial efficacies can dissipate their positive influence on the evolution of the infection. At the level of the evolution of the Gini index, in the presence of reinfection, it appears even more evident that inequalities appear for large times in the presence of mild vaccinations. Nevertheless, in transient regimes, the higher possibility of investing wealth for vaccinated agents may create temporary inequalities.
