# From agent-based models to the macroscopic description of fake-news spread:
the role of competence in data-driven applications

J. Franceschi<EMAIL_ADDRESS>Department of
Mathematics “F. Casorati”, University of Pavia, Italy  L. Pareschi
<EMAIL_ADDRESS>Department of Mathematics and Computer Science,
University of Ferrara, Italy  M. Zanella<EMAIL_ADDRESS>Department of
Mathematics “F. Casorati”, University of Pavia, Italy

###### Abstract

Fake news spreading, with the aim of manipulating individuals’ perceptions of
facts, is now recognized as a major problem in many democratic societies. Yet,
to date, little has been understood about how fake news spreads on social
networks, what the influence of the education level of individuals is, when
fake news is effective in influencing public opinion, and what interventions
might be successful in mitigating their effect. In this paper, starting from
the recently introduced kinetic multi-agent model with competence by the first
two authors, we propose to derive reduced-order models through the notion of
social closure in the mean-field approximation that has its roots in the
classical hydrodynamic closure of kinetic theory. This approach allows to
obtain simplified models in which the competence and learning of the agents
maintain their role in the dynamics and, at the same time, the structure of
such models is more suitable to be interfaced with data-driven applications.
Examples of different Twitter-based test cases are described and discussed.

Keywords: fake news spreading, learning dynamics, agent-based models, kinetic
models, social closure, competence, data uncertainty

###### Contents

  1. 1 Introduction
  2. 2 Kinetic models, competence and fake news spreading
    1. 2.1 Competence and learning in multi-agent systems
    2. 2.2 Fake news spreading in presence of a social feature
    3. 2.3 Asymptotic states of the learning process
  3. 3 Reduced order models for fake news spread with competence
    1. 3.1 Social closure with a strong competence-based contact function
    2. 3.2 Social closure with a weak competence-based contact function
  4. 4 Examples and applications
    1. 4.1 Test 1: Validation of the social closure
    2. 4.2 Test 2: A data driven application to Twitter
      1. 4.2.1 Test 2A: Fitting the model to data
      2. 4.2.2 Test 2B: Forecasting under data uncertainties
    3. 4.3 Test 3: Competence background in misinformation
  5. Acknowledgments

##  1 Introduction

Since the 2016 U.S. presidential election, and more recently the COVID-19
infodemic, fake news on social networks, intended to manipulate users’
perceptions of events, has been recognized as a fundamental problem in open
societies. As fake news proliferate, disinformation threatens democracy and
efficient governance. In particular, there is empirical evidence that fake
news spreads significantly “faster, deeper, and more widely” than real news
[37]. In the same study, it is also highlighted that the phenomenon is not due
to robotic automatisms of news dissemination but to the actions of human
beings sharing the news without the ability to identify misinformation.

It is therefore of fundamental importance the construction of mathematical
models capable of describing such scenarios and with a structure simple enough
to be interfaced with data available, for example from social networks, but
still embedding the specific features related to the ability of individuals in
detecting the piece of false information.

In recent years, compartmental models inspired by epidemiology have been used
fruitfully to study spreading phenomena of rumors and hoaxes. For instance,
following the pioneering work of Daley & Kendall [11], in [23] SIR-type models
are used in conjunction with dynamical trust rates that account for the
different spreading rates in a network. Those traditional models were
elaborated in [7], where the authors consider also the impact of online groups
in feeding the rumor growth once it has started.

Alongside these approaches there are more data-driven works. In this field,
Twitter has been gaining consensus as a powerful source of useful and
structured information. A recent example in this direction can be found in
[27], that focuses on fake news dissemination on the platform using a two-
phase model, where fake news initially spread as novel news story and after a
correction time they are paired with a competitive narrative which describes
the news as fake in the first place.

Twitter data in conjunction with epidemiological models have already been used
to study the spread of rumors and fake news by several authors [16, 17, 26,
10], where SIS and SEIZ compartmental models were employed to fit the data of
the evolution of different news. Mounting experimental evidence highlights the
strong link between digital media literacy and possibility to reliably
identify the quality of online information. This connection has been early
identified by communication scientists [20] and later confirmed by
experimental studies, see e.g. [22, 24].

In [19], starting from an agent-based model for the dissemination of fake news
in presence of competence, using the tools of kinetic theory, in the limit of
a large number of agents, novel mathematical models were proposed and
discussed. Previously, kinetic models that include the role of competence or
knowledge had been proposed in [5, 29, 31]. The behavior of a social system
composed by a large number of interacting agents has been studied in the case
of opinion formation [3, 9, 14, 15, 34] and more recently epidemiological
dynamics [1, 2, 12]. We refer to [28] for an introduction to the subject.

The compartmental structure of the model for fake-news spreading in presence
of competence introduced in [19] is composed by four groups of individuals:
the susceptible (S) agents—defined as the ones who are unaware of the fake
news; the exposed (E) agents—those who know the news but still have not decide
whether to spread it or not; the infectious (I) agents—who actively divulge
and finally the skeptical or removed (R) agents—those who are aware of the
news but choose to not spread it. On a population divided among such
categories, there is also a social structure based on an additional time
evolving variable that measures the _competence_ level of the agents. Although
the model has shown the capacity to correctly describe the role of competence
in the dynamics of fake-news, its mathematical structure based on kinetic
partial differential equations is generally too complex to be interfaced with
the available data.

In an attempt to address this problem, in the present work by exploiting the
knowledge of the equilibrium states of the corresponding mean-field model we
derived reduced order macroscopic models based on ordinary differential
equations in which, however, the role of competence continues to be present.
The new social models, thanks to their simpler structure, are more suitable
for data-driven applications. We emphasize that the methodology here adopted
is quite general and that in principle points the way to introducing
additional social characteristics of individuals into tractable mathematical
models in terms of structural complexity.

The rest of the manuscript is organized according to the following sections.
In Section 2, we recall the basic concepts of the kinetic model for describing
the spread of fake-news in the presence of competence. Next, in Section 3,
using the local equilibrium states of the competence we derive reduced order
models that depend on the specific shape of the interaction function. Section
4 is devoted to presenting a series of numerical experiments in which we first
validate the model and then consider data-driven applications based on
Twitter. In the last section, a series of final considerations are reported.

##  2 Kinetic models, competence and fake news spreading

In this section we present a model for the description of the spreading of
fake news in a society characterized by a heterogeneous competence of agents.
Our starting point is the compartmental kinetic approach recently proposed in
[19]. We suppose that the system of agents can be divided in the following
epidemiologically relevant states: susceptible (S) agents are the ones that
are unaware of fake news, we further denote as exposed (E) the agents that
encountered the fake news but have still to spread them, infectious (I) agents
are the real spreader and, finally, the removed (R) agents are not actively
engaged in the spread of misinformation. In the following we indicate with
$\mathcal{C}=\\{S,E,I,R\\}$ the set of epidemiological compartments.

Aiming to incorporate the effects of personal competence on the fake news
dynamics, we stick to a simple mathematical setting where the state of the
individuals in each compartment, at any time $t\geq 0$, is characterized by
the sole competence level $x\in\mathbb{R}_{+}$. Hence, we denote by
$f_{S}=f_{S}(x,t)$, $f_{E}=f_{E}(x,t)$, $f_{I}=f_{I}(x,t)$, and
$f_{R}=f_{R}(x,t)$ the distribution of competence at time $t\geq 0$ of
susceptible, exposed, infectious and removed individuals, respectively. We
neglect natality and mortality dynamics since we can consider a short time
dynamic where nobody enters or leaves it during the spreading of the fake
news. This assumption can be justified based on the average lifespan of fake
news. Therefore, we can fix the total distribution of competence of a society
to be a probability density for all $t\geq 0$

$\int_{\mathbb{R}_{+}}\left(f_{S}(x,t)+f_{E}(x,t)+f_{I}(x,t)+f_{R}(x,t)\right)\,\mathrm{d}x=1,\quad
t>0,$

Consequently, the quantities

$\displaystyle S(t)$
$\displaystyle=\int_{\mathbb{R}_{+}}f_{S}(x,t)\,\mathrm{d}x,$ $\displaystyle
E(t)$ $\displaystyle=\int_{\mathbb{R}_{+}}f_{E}(x,t)\,\mathrm{d}x,$
$\displaystyle I(t)$
$\displaystyle=\int_{\mathbb{R}_{+}}f_{I}(x,t)\,\mathrm{d}x,$ $\displaystyle
R(t)$ $\displaystyle=\int_{\mathbb{R}_{+}}f_{R}(x,t)\,\mathrm{d}x$

denote the fractions of the population that are susceptible, exposed,
infected, or recovered respectively at time $t\geq 0$. We also denote with
$m_{J}^{p}(t)$ the moment of the distribution $f_{J}(x,t)$, $J\in\mathcal{C}$,
of order $p\geq 0$

$\displaystyle m^{p}_{S}(t)$
$\displaystyle=\frac{1}{S(t)}\int_{\mathbb{R}_{+}}x^{p}f_{S}(x,t)\,\mathrm{d}x,$
$\displaystyle m^{p}_{E}(t)$
$\displaystyle=\frac{1}{E(t)}\int_{\mathbb{R}_{+}}x^{p}f_{E}(x,t)\,\mathrm{d}x,$
$\displaystyle m^{p}_{I}(t)$
$\displaystyle=\frac{1}{I(t)}\int_{\mathbb{R}_{+}}x^{p}f_{I}(x,t)\,\mathrm{d}x,$
$\displaystyle m^{p}_{R}(t)$
$\displaystyle=\frac{1}{R(t)}\int_{\mathbb{R}_{+}}x^{p}f_{R}(x,t)\,\mathrm{d}x.$

Unambiguously we will indicate with $m_{J}(t)$, $J\in\mathcal{C}$, the mean
values corresponding to $p=1$.

###  2.1 Competence and learning in multi-agent systems

Drawing inspiration from seminal models for multi-agent systems in presence of
personal competence [29, 31] we introduce a binary interaction term expressing
two different processes:

  1. it)

learning processes by less competent agents that can learn from the more
competent ones

  2. it)

the competence evolution depends by a social background in which individuals
grow.

The dynamics described at point 1 can be easily sketched by the following
process: if two agents belonging to compartment $H,J\in\mathcal{C}$ and
characterized by competence levels $x,x_{*}\in\mathbb{R}_{+}$ meet, their
post-interaction competence is given by

$\displaystyle\left\\{\begin{aligned}
x^{\prime}&=(1-\lambda_{H}(x))x+\lambda_{CJ}(x)x_{\ast}+\eta_{HJ}x\\\
x_{\ast}^{\prime}&=(1-\lambda_{J}(x_{\ast}))x_{\ast}+\lambda_{CH}(x_{\ast})x+\eta_{JH}x_{\ast},\end{aligned}\right.$
$\displaystyle H,J\in\mathcal{C}$ (1)

where $\lambda_{H}(\cdot)$, $H\in\mathcal{C}$, quantify the amount of
competence lost by individuals of compartment $H$ by the natural process of
forgetfulness and the parameter $\lambda_{CH}$, $H\in\mathcal{C}$, models the
competence gained through the interaction with members of the class $J$, with
$J\in\mathcal{C}$. A possible choice for $\lambda_{CJ}(x)$ is
$\lambda_{CJ}(x)=\lambda_{CJ}\chi(x\geq\bar{x})$, where $\chi(\cdot)$ is the
characteristic function and $\bar{x}\in X$ a minimum level of competence
required to the agents for increasing their own skills by interactions. In (1)
$\eta_{HJ}$ and $\eta_{JH}$ are centered iid random variable such that,
denoting by $\left\langle\cdot\right\rangle$ their expectation, we have
$\left\langle\eta_{HJ}^{2}\right\rangle=\left\langle\eta_{JH}^{2}\right\rangle=\sigma^{2}_{HJ}$.

We suppose that the process defined in $(ii)$ takes place in a different time
scale from the one of interactions between agents. In particular, unlike [19]
we assume that the time scale of online interactions for competence formation
is faster than interactions with the social background. To this end, we will
consider advection terms that will be defined in the next section.

######  Remark 1.

It is reasonable to assume that both the processes of gain and loss of
competence from the interaction with other agents in (1) are bounded by zero.
Therefore we suppose that if $J,H\in\\{S,E,I,R\\}$, and if
$\lambda_{J}\in[\lambda_{J}^{-},\lambda_{J}^{+}]$, with $\lambda_{J}^{-}>0$
and $\lambda_{J}^{+}<1$, and $\lambda_{CJ}(x)\in[0,1]$ then $\eta_{HJ}$ may,
for example, be uniformly distributed in
$[-1+\lambda_{J}^{+},1-\lambda_{J}^{+}]$.

###  2.2 Fake news spreading in presence of a social feature

Following [19] we choose to describe the dissemination of fake news through a
population of agents via a kinetic compartmental model. In this setting the
description of the sole spreading dynamics can be illustrated by the following
system of ODEs

$\displaystyle\left\\{\begin{aligned} \frac{dS}{dt}&=-\beta
SI+(1-\alpha)\gamma I\\\ \frac{dE}{dt}&=\beta SI-\delta E\\\
\frac{dI}{dt}&=(1-\eta)\delta E-\gamma I\\\ \frac{dR}{dt}&=\eta\delta
E+\alpha\gamma I\end{aligned}\right.$ (2) $\displaystyle S+E+I+R=1.$

Borrowing from the consolidated epidemiological tradition, we will refer to it
as the SEIR model. System (2) describes the evolution of the mass fractions of
the population that belongs to each compartment $J\in\mathcal{C}$ for each
time $t\geq 0$. The parameters appearing in system (2) are presented in Table
1. Also a schematic representation of system (2) is given in Figure 2.2. The
last equation of system (2) translates the fact that—as specified at the
beginning of the Section—the total mass of the population is preserved.

The combination of the learning mechanisms presented in the previous
subsection together with the spreading of the fake news is described by the
following kinetic model:

$\left\\{\begin{aligned} \frac{\partial f_{S}(x,t)}{\partial
t}&=-K(f_{S},f_{I})(x,t)+(1-\alpha(x))\gamma(x)f_{I}(x,t)+\dfrac{1}{\varepsilon}\sum_{\mathclap{J\in\mathcal{C}}}\mathcal{Q}_{SJ}(f_{S},f_{J})(x,t)\\\
&\phantom{=}+\frac{\partial}{\partial
x}\left[\left(\lambda_{S}x-\lambda_{BS}m_{B}\right)f_{S}(x,t)\right],\\\
\frac{\partial f_{E}(x,t)}{\partial
t}&=K(f_{S},f_{I})(x,t)-\delta(x)f_{E}(x,t)+\dfrac{1}{\varepsilon}\sum_{\mathclap{J\in\mathcal{C}}}\mathcal{Q}_{EJ}(f_{E},f_{J})(x,t)\\\
&\phantom{=}+\frac{\partial}{\partial
x}\left[\left(\lambda_{E}x-\lambda_{BE}m_{B}\right)f_{E}(x,t)\right],\\\
\frac{\partial f_{I}(x,t)}{\partial
t}&=\delta(x)(1-\eta(x))f_{E}(x,t)-\gamma(x)f_{I}(x,t)+\dfrac{1}{\varepsilon}\sum_{\mathclap{J\in\mathcal{C}}}\mathcal{Q}_{\mathit{IJ}}(f_{I},f_{J})(x,t)\\\
&\phantom{=}+\frac{\partial}{\partial
x}\left[\left(\lambda_{I}x-\lambda_{BI}m_{B}\right)f_{I}(x,t)\right],\\\
\frac{\partial f_{R}(x,t)}{\partial
t}&=\delta(x)\eta(x)f_{E}(x,t)+\alpha(x)\gamma(x)f_{I}(x,t)+\dfrac{1}{\varepsilon}\sum_{\mathclap{J\in\mathcal{C}}}\mathcal{Q}_{RJ}(f_{R},f_{J})(x,t)\\\
&\phantom{=}+\frac{\partial}{\partial
x}\left[\left(\lambda_{R}x-\lambda_{BR}m_{B}\right)f_{R}(x,t)\right],\end{aligned}\right.$
(3)

where the parameter $\varepsilon$ describes the intensity of the interactions.

Parameter | Definition
---|---
$\beta$ | contact rate between susceptible and infected individuals
$1/\delta$ | average decision time on whether or not to spread fake news
$\eta$ | probability of deciding not to spread fake news
$1/\gamma$ | average duration of a fake news
$\alpha$ | probability of remembering fake news
Table 1: Parameters definition in the SEIR model (3).

In (3) the functional

$K(f_{S},f_{I})(x,t)=f_{S}(x,t)\int_{\mathbb{R}^{+}}\kappa(x,x_{\ast})f_{I}(x_{\ast},t)\,\mathrm{d}x_{\ast}$
(4)

is the local incidence rate and $\kappa(x,x_{\ast})$ is a nonnegative contact
function measuring the impact of competence in the spreading of fake news.
This function is decreasing with respect to the competences $x,x_{*}\geq 0$ of
the population of susceptible and infected agents. In the following we will
investigate the macroscopic effects of the following two choices of
$\kappa(x,x_{\ast})$

  1. it)

Strong competence-based contact function
$\kappa(x,x_{\ast})=\beta/(x\,x_{*})$, with $\beta>0$,

  2. it)

Weak competence-based contact function $\kappa(x,x_{\ast})=\beta
e^{-x-x_{*}}$, $\beta>0$.

The two functions are both decreasing but have strong differences for
$x,x_{*}\ll 1$. Indeed, since 1 is not limited for small competences it
enforces the spreading of fake news among less competent agents compared with
$(B)$. Indeed, the function in 2 is bounded in $\mathbb{R}_{+}$. We further
remark that individuals have the highest rates of contact with people
belonging to the same social class, and thus with a similar level of
competence.

######  Remark 2.

In (3) the parameters are considered to be dependent on the competence level
$x$, in general. This is to reflect the fact that competence plays a role in
the dissemination of fake news.

SERI$\beta$$(1-\eta)\delta$$\eta\delta$

$\alpha\gamma$

$(1-\alpha)\gamma$

Figure 1: Dissemination dynamics for the SEIR model (2)

Furthermore, the operators $\mathcal{Q}_{HJ}(f_{H},f_{J})(x,t)$,
$J\in\mathcal{C}$, describe the binary collisions (1) and they determine the
thermalization of the distribution of competence characterizing the $J$th
compartment. The advection terms in (3) come models the influence of the
social background on the competence dynamics. It is worth to observe that the
evolution of mass fractions $J(t)$ obeys the classical SEIR model with
reinfection (2) by choosing $\mathcal{Q}_{HJ}\equiv 0$ and
$\kappa(x,x_{*})=\beta>0$. This would correspond in considering the spreading
a fake news independent of the competence level of a system of agents.

In more details, we will consider the operators $\mathcal{Q}_{HJ}$ as integral
operators that modify the competence distribution through repeated
interactions of type (1) among individuals. We can fruitfully define the
introduced operators in weak form as follows

$\int_{\mathbb{R}^{+}}\varphi(x)\mathcal{Q}_{HJ}(f_{H},f_{J})\,\mathrm{d}x=\left\langle\int_{\mathbb{R}^{2}_{+}}\bigl{(}\varphi(x^{\prime})-\varphi(x)\bigr{)}\,f_{H}(x,t)f_{J}(x_{\ast},t)\,\mathrm{d}x_{\ast}\mathrm{d}x\right\rangle,$
(5)

where $\varphi(\cdot)$ is a test function and where the brackets
$\langle\cdot\rangle$ indicate the expectation with respect to the random
variables $\eta_{HJ},\tilde{\eta}_{HJ}$.

In the model (3) the function $\gamma(x)>0$ determines the duration of the
fake news and can be strongly influenced by the competence level of the
spreader. Furthermore, the function $\delta(x)>0$ is related to the average
time that an agent eventually spend before the diffusion of a fake news such
that people with high competence invest more time in checking information
reliability, and $\eta(x)\in[0,1]$ characterizes individuals’ decision to
spread fake news. The function $\alpha(x)\in[0,1]$ describes the probability
to remember fake news and can be thought less influenced by the competence
variable. In Table 1 we summarize all the introduced parameters.

###  2.3 Asymptotic states of the learning process

We focus now on the learning dynamics introduced in model (3) whose evolution
is given by the nonlinear operators $\mathcal{Q}_{HJ}(f_{H},f_{J})$,
$H,J\in\mathcal{C}$, defined in (5). We concentrate in particular on the
analysis of asymptotic states of the learning dynamics undergoing elementary
interactions (1). We are therefore interested in the asymptotic distribution
of the Boltzmann-type model

$\begin{split}\dfrac{d}{dt}\int_{\mathbb{R}_{+}}\varphi(x)f_{H}(x,t)dx&=\sum_{J\in\mathcal{C}}\int_{\mathbb{R}^{2}_{+}}\left\langle\varphi(x^{\prime})-\varphi(x)\right\rangle\,f_{H}(x,t)f_{J}(x_{\ast},t)\,\mathrm{d}x_{\ast}\mathrm{d}x\\\
&=\dfrac{1}{2}\sum_{J\in\mathcal{C}}\int_{\mathbb{R}^{2}_{+}}\left\langle\varphi(x^{\prime})+\varphi(x_{*}^{\prime})-\varphi(x)-\varphi(x_{*})\right\rangle\,f_{H}(x,t)f_{J}(x_{\ast},t)\,\mathrm{d}x_{\ast}\mathrm{d}x.\\\
\end{split}$ (6)

It is easily observed that if $\varphi(x)=1$ the mass is conserved in (6)
corresponding to the conservation of the total number of agents. If
$\varphi(x)=x$ in (6) we obtain the evolution of the average competence in
each compartment that is not conserved in time

$\begin{split}\dfrac{d}{dt}(H(t)m_{H}(t))&=\sum_{J\in\mathcal{C}}\int_{\mathbb{R}_{+}^{2}}\left\langle
x^{\prime}-x\right\rangle
f_{H}(x,t)f_{J}(x_{*},t)\mathrm{d}x\,\mathrm{d}x_{*}\\\
&=H(t)\sum_{J\in\mathcal{C}}J(t)(\lambda_{CJ}m_{J}(t)-\lambda_{H}m_{H}(t)),\end{split}$

and the total competence is conserved

$\dfrac{d}{dt}\sum_{H\in\mathcal{C}}\int_{\mathbb{R}_{+}}xf_{H}(x,t)\mathrm{d}x=0.$

Since the steady state solution of (6) is difficult to obtain, we can formally
derive a simplified Fokker-Planck model in which the study of the asymptotic
properties is much easier. To this end, we introduce the following quasi-
invariant scaling of the relevant parameter of the binary scheme (1) given by

$\lambda_{H}\rightarrow\tau\lambda_{H},\qquad\lambda_{CH}\rightarrow\tau\lambda_{CH},\qquad\sigma_{HJ}^{2}\rightarrow\tau\sigma_{HJ}^{2},$
(7)

with $\tau>0$. It is worth to mention that the introduced scaling is inspired
by the so-called grazing collision limit of the Boltzmann equation, see [6,
36]. In the context of multi-agent systems this scaling has been introduced in
[8, 33].

In the introduced regime of parameters the interactions become quasi-
invariant, in the sense that the post-interaction competences
$(x^{\prime},x_{*}^{\prime})$ are such that $x^{\prime}-x$ and
$x_{*}^{\prime}-x_{*}$ are small for $\tau\ll 1$. Hence, assuming
$\varphi\in\mathcal{C}_{0}$, we can perform the following Taylor expansion

$\begin{split}\varphi(x^{\prime})-\varphi(x)&=(x^{\prime}-x)\dfrac{\mathrm{d}}{\mathrm{d}x}\varphi(x)+\dfrac{1}{2}(x^{\prime}-x)^{2}\dfrac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\varphi(x)+\dfrac{1}{6}(x^{\prime}-x)^{3}\dfrac{\mathrm{d}^{3}}{\mathrm{d}x^{2}}\varphi(\bar{x}),\end{split}$

with $\bar{x}\in(\min\\{x,x^{\prime}\\},\max\\{x,x^{\prime}\\})$. Hence, in
the time scale $t/\tau$ we have

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathbb{R}_{+}}\varphi(x)f_{H}(x,t)\mathrm{d}x\\\
=\sum_{J\in\mathcal{C}}\int_{\mathbb{R}^{2}_{+}}(-\lambda_{H}(x)x+\lambda_{CJ}(x)x_{*})\frac{\mathrm{d}\varphi(x)}{\mathrm{d}x}f_{H}(x,t)f_{J}(x_{*},t)\,\mathrm{d}x_{*}\,\mathrm{d}x\\\
\quad+\sum_{J\in\mathcal{C}}\frac{\sigma_{HJ}^{2}}{2}\int_{\mathbb{R}_{+}^{2}}x^{2}\frac{\mathrm{d}^{2}\varphi(x)}{\mathrm{d}x^{2}}f_{H}(x,t)f_{J}(x_{*},t)\,\mathrm{d}x_{*}\,\mathrm{d}x+\sum_{J\in\mathcal{C}}R_{\varphi}(f_{H},f_{J}),$

where we exploited the fact that $\left\langle\eta_{HJ}\right\rangle=0$ and we
have defined the sum of reminder terms

$\begin{split}R_{\tau}(f_{H},f_{J})&=\frac{1}{2}\int_{\mathbb{R}_{+}^{2}}\tau(-\lambda_{H}(x)x+\lambda_{CJ}(x)x_{*})^{2}\frac{\mathrm{d}^{2}\varphi(x)}{\mathrm{d}x^{2}}f_{H}(x,t)f_{J}(x_{*},t)\mathrm{d}x\,\mathrm{d}x_{*}\\\
&\phantom{{}={}}+\frac{1}{6}\int_{\mathbb{R}_{+}^{2}}\frac{\left\langle(-\lambda_{H}(x)x+\lambda_{CJ}(x)x_{*}+\eta_{HJ}x)^{3}\right\rangle}{\tau}\frac{\mathrm{d}^{3}\varphi(x)}{\mathrm{d}x^{3}}f_{H}(x,t)f_{J}(x_{*},t)\mathrm{d}x\,\mathrm{d}x_{*}\end{split}$

We may observe that, assuming
$\left\langle|\eta_{HJ}|^{3}\right\rangle<+\infty$, then we may write
$\eta_{HJ}=\sqrt{\sigma^{2}}\tilde{\eta}_{HJ}$, where we introduced the
centered random variable $\tilde{\eta}_{HJ}$ with unitary variance and such
that $\left\langle|\tilde{\eta}|^{3}\right\rangle<+\infty$. Therefore,
$\left\langle|\eta_{HJ}|^{3}\right\rangle=(\sigma^{2})^{3/2}\left\langle|\tilde{\eta}_{HJ}|^{3}\right\rangle$
and, under the scaling (7), we get
$\left\langle|\eta_{HJ}|^{3}\right\rangle\sim\tau^{3/2}\sigma^{3/2}$ . Hence,
under the above assumption, proceeding as in [8] we can prove that for
$\tau\rightarrow 0^{+}$

$|R_{\varphi}(f_{H},f_{J})(x,t)|\rightarrow 0.$

Therefore in the new time scale, for $\tau\rightarrow 0^{+}$ and under the
quasi-invariant scaling (7), we can show that the solution of model (6)
converges to

$\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathbb{R}_{+}}\varphi(x)f_{H}(x,t)\mathrm{d}x\\\
=\sum_{J\in\mathcal{C}}\int_{\mathbb{R}^{2}_{+}}(-\lambda_{H}(x)x+\lambda_{CJ}(x)x_{*})\frac{\mathrm{d}\varphi(x)}{\mathrm{d}x}f_{H}(x,t)f_{J}(x_{*},t)\,\mathrm{d}x_{*}\,\mathrm{d}x\\\
\quad+\sum_{J\in\mathcal{C}}\dfrac{\sigma_{HJ}^{2}}{2}\int_{\mathbb{R}_{+}^{2}}x^{2}\frac{\mathrm{d}^{2}\varphi(x)}{\mathrm{d}x^{2}}f_{H}(x,t)f_{J}(x_{*},t)\,\mathrm{d}x_{*}\,\mathrm{d}x.$
(8)

Integrating back by parts we have obtained

$\partial_{t}f_{H}(x,t)=\partial_{x}\left[(\lambda_{H}(x)x-\sum_{J\in\mathcal{C}}\lambda_{CJ}(x)J(t)m_{J}(t))f_{H}(x,t)+\dfrac{\sigma^{2}}{2}\partial_{x}(x^{2}f_{H}(x,t))\right],$
(9)

with $\sum_{J\in\mathcal{C}}\sigma_{HJ}^{2}=\sigma^{2}$, coupled with the
following boundary conditions

$(\lambda_{H}(x)x-\sum_{J\in\mathcal{C}}\lambda_{CJ}(x)J(t)m_{J}(t))f_{H}(x,t)+\dfrac{\sigma^{2}}{2}\partial_{x}(x^{2}f_{H}(x,t))\Big{|}_{x=0}=0,$

and

$\dfrac{\sigma^{2}}{2}\partial_{x}(x^{2}f_{H}(x,t))\Big{|}_{x=0}=0,$

for all $H\in\mathcal{C}$. Assuming then $\lambda_{CJ}=\lambda_{H}=\lambda$
independent by $x\in\mathbb{R}_{+}$ for all $J,H\in\mathcal{C}$ the steady
states $f_{H}^{\infty}(x)$, $H\in\mathcal{C}$ are solution of

$\lambda(x-m)f_{H}^{\infty}(x)+\dfrac{\sigma^{2}}{2}\partial_{x}(x^{2}f_{H}^{\infty}(x))=0,$

where $m=\sum_{J\in\mathcal{C}}J(t)m_{J}(t)$ is a conserved quantity as we
already observed. Hence, we obtain that the large time distribution is an
inverse Gamma

$f_{H}^{\infty}(x)=H\dfrac{k^{\mu}}{\Gamma(\mu)}\dfrac{e^{-k/x}}{x^{1+\mu}},$

where

$\mu=1+\dfrac{2\lambda}{\sigma^{2}},\qquad k=(\mu-1)m.$

Now, we highlight that for $t\rightarrow+\infty$ we have
$f_{E}(x,t),f_{I}(x,t)\rightarrow 0$, which means that

$f^{\infty}(x)=f_{S}^{\infty}(x)+f_{R}^{\infty}(x).$

In view of $S^{\infty}+R^{\infty}=1$ we conclude that under the introduced
assumptions

$f_{S}^{\infty}(x)=S^{\infty}f^{\infty},\qquad
f_{R}^{\infty}=(1-S^{\infty})f^{\infty}.$

##  3 Reduced order models for fake news spread with competence

Once we have characterized the equilibrium distribution of the transition
operators $\mathcal{Q}_{HJ}(\cdot,\cdot)$, with $H,J\in\mathcal{C}$, we can
study the complete system (3). The aim of this section is the definition of
observable macroscopic equations of the introduced kinetic model.

Integrating both sides of (3) with respect to $x\in\mathbb{R}_{+}$ and
recalling that the introduced operators are mass and momentum preserving, we
obtain the following system for the evolution of the mass fractions $J(t)$,
$J\in\mathcal{C}$

$\left\\{\begin{aligned}
\frac{\mathrm{d}S(t)}{\mathrm{d}t}&=-\int_{\mathbb{R}^{2}_{+}}\kappa(x,x_{*})f_{S}(x,t)f_{I}(x_{*},t)dx\,dx_{*}+(1-\alpha)\gamma
I(t),\\\
\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\int_{\mathbb{R}^{2}_{+}}\kappa(x,x_{*})f_{S}(x,t)f_{I}(x_{*},t)dx\,dx_{*}-\delta
E(t),\\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}&=\delta(1-\eta)E(t)-\gamma
I(t),\\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}&=\delta\eta E(t)+\alpha\gamma
I(t),\end{aligned}\right.$ (10)

whereas for the momentum we get

$\left\\{\begin{aligned}
\frac{\mathrm{d}(S(t)m_{S}(t))}{\mathrm{d}t}&=-\int_{\mathbb{R}^{2}_{+}}x\kappa(x,x_{*})f_{S}(x,t)f_{I}(x_{*},t)dx\,dx_{*}+(1-\alpha)\gamma
I(t)m_{I}(t),\\\
\frac{\mathrm{d}(E(t)m_{E}(t))}{\mathrm{d}t}&=\int_{\mathbb{R}^{2}_{+}}x\kappa(x,x_{*})f_{S}(x,t)f_{I}(x_{*},t)dx\,dx_{*}-\delta
E(t),\\\
\frac{\mathrm{d}(I(t)m_{I}(t))}{\mathrm{d}t}&=\delta(1-\eta)E(t)m_{E}(t)-\gamma
I(t)m_{I}(t),\\\ \frac{\mathrm{d}(R(t)m_{R}(t))}{\mathrm{d}t}&=\delta\eta
E(t)m_{E}(t)+\alpha\gamma I(t)m_{I}(t).\end{aligned}\right.$ (11)

We can observe that the obtained system is not closed since the evolution of
mass fractions $J(t)$ and of the momentum depend on the evolution of the
distribution functions $f_{J}(x,t)$. The closure of the obtained system can be
obtained by formally resorting to a limit procedure. Indeed, assuming that the
time scale involved in the process of competence formation is $\varepsilon\ll
1$, we have a fast learning process of the system of agents with respect to
the evolution of the spreading of fake news. Therefore, for $\varepsilon\ll 1$
the distribution function $f_{J}(x,t)$ reaches fast the inverse Gamma
equilibrium with mass fractions $J(t)$ and local mean values $m_{J}(t)$.

In the following we obtain two different set of macroscopic equations in
relation with the considered contact rate function $\kappa(x,x_{*})$.

###  3.1 Social closure with a strong competence-based contact function

We consider the case $(A)$ introduced in Section 2.2 corresponding to a strong
competence-based contact function defined by
$\kappa(x,x_{*})=\frac{\beta}{xx_{*}}$, $\beta>0$. We have

$\left\\{\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}&=-\beta
H_{S}(t)S(t)H_{I}(t)I(t)+(1-\alpha)\gamma I(t),\\\
\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\beta H_{S}(t)S(t)H_{I}(t)I(t)-\delta
E(t),\\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}&=\delta(1-\eta)E(t)-\gamma
I(t),\\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}&=\delta\eta E(t)+\alpha\gamma
I(t),\end{aligned}\right.$ (12)

where

$H_{J}(t)=\int_{\mathbb{R}^{+}}\frac{1}{x}f_{J}(x,t)\,\mathrm{d}x.$ (13)

Therefore, in the limit $\varepsilon\rightarrow 0^{+}$ we can plug
$f_{J}^{\infty}(x)$ in (13) which becomes

$H_{J}(t)=\frac{\mu}{\mu-1}\frac{1}{m_{J}(t)},$ (14)

thanks to the properties of the inverse Gamma distribution, leading to

$\left\\{\begin{aligned}
\frac{\mathrm{d}S(t)}{\mathrm{d}t}&=-\beta\biggl{(}\frac{\mu}{\mu-1}\biggr{)}^{2}\frac{S(t)I(t)}{m_{S}(t)m_{I}(t)}+(1-\alpha)\gamma
I(t),\\\
\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\beta\biggl{(}\frac{\mu}{\mu-1}\biggr{)}^{2}\frac{S(t)I(t)}{m_{S}(t)m_{I}(t)}-\delta
E(t),\\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}&=\delta(1-\eta)E(t)-\gamma
I(t),\\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}&=\delta\eta E(t)+\alpha\gamma
I(t).\end{aligned}\right.$ (15)

Next, looking at (11), recalling that under the hypothesis that
$\lambda_{J}=\lambda_{CJ}$ for $J\in\mathcal{C}$, the knowledge exchange
operator also preserves momentum, we have the following system of equations

$\left\\{\begin{aligned} \frac{\mathrm{d}[m_{S}(t)S(t)]}{\mathrm{d}t}&=-\beta
S(t)\frac{\mu}{\mu-1}\frac{I(t)}{m_{I}(t)}+(1-\alpha)\gamma
m_{I}(t)I(t)+\lambda(m_{B}-m_{S}(t))S(t),\\\
\frac{\mathrm{d}[m_{E}(t)E(t)]}{\mathrm{d}t}&=\beta
S(t)\frac{\mu}{\mu-1}\frac{I(t)}{m_{I}(t)}-\delta
m_{E}(t)E(t)+\lambda(m_{B}-m_{E}(t))E(t),\\\
\frac{\mathrm{d}[m_{I}(t)I(t)]}{\mathrm{d}t}&=\delta(1-\eta)m_{E}(t)E(t)-\gamma
m_{I}(t)I(t)+\lambda(m_{B}-m_{I}(t))I(t),\\\
\frac{\mathrm{d}[m_{R}(t)R(t)]}{\mathrm{d}t}&=\delta\eta
m_{E}(t)E(t)+\alpha\gamma
m_{I}(t)I(t)+\lambda(m_{B}-m_{R}(t))R(t),\end{aligned}\right.$ (16)

which, using the fact that

$\frac{\mathrm{d}[m_{J}(t)J(t)]}{\mathrm{d}t}=J(t)\frac{\mathrm{d}m_{J}(t)}{\mathrm{d}t}+m_{J}(t)\frac{\mathrm{d}J(t)}{\mathrm{d}t},$

implies

$\left\\{\begin{aligned}
\frac{\mathrm{d}m_{S}(t)}{\mathrm{d}t}&=\beta\biggl{(}\frac{\mu}{\mu-1}\biggr{)}\frac{I(t)}{m_{I}(t)}\frac{1}{\mu-1}+(1-\alpha)\gamma\frac{I(t)}{S(t)}\biggl{[}m_{I}(t)-m_{S}(t)\biggr{]}\\\
&\phantom{{}={}}+\lambda(m_{B}-m_{S}(t)),\\\
\frac{\mathrm{d}m_{E}(t)}{\mathrm{d}t}&=\beta\biggl{(}\frac{\mu}{\mu-1}\biggr{)}\frac{S(t)I(t)}{E(t)m_{I}(t)}\biggl{[}1-\frac{\mu}{\mu-1}\frac{m_{E}(t)}{m_{S}(t)}\biggr{]}+\lambda(m_{B}-m_{E}(t)),\\\
\frac{\mathrm{d}m_{I}(t)}{\mathrm{d}t}&=\delta(1-\eta)\frac{E(t)}{I(t)}(m_{E}(t)-m_{I}(t))+\lambda(m_{B}-m_{I}(t)),\\\
\frac{\mathrm{d}m_{R}(t)}{\mathrm{d}t}&=\delta\eta\frac{E(t)}{R(t)}(m_{E}(t)-m_{R}(t))+\alpha\gamma\frac{I(t)}{R(t)}(m_{I}(t)-m_{R}(t))\\\
&\phantom{{}={}}+\lambda(m_{B}-m_{R}(t)),\end{aligned}\right.$ (17)

that is, we obtained a closed system of eight ordinary differential equations
(15), (17).

###  3.2 Social closure with a weak competence-based contact function

If, instead, we consider the case 2 of Section 2.2, corresponding to the weak
competence-based contact function defined by
$\kappa(x,y)=\mathrm{e}^{-x}\mathrm{e}^{-y}$, it is possible to write

$\tilde{H}_{J}(t)=\int_{\mathbb{R}_{+}}e^{-x}f_{J}(x,t)dx.$ (18)

As discussed in Section 3.1, in the limit $\varepsilon\rightarrow 0^{+}$ we
may plug the asymptotic distribution $f_{J}^{\infty}$ of the Fokker-Planck
model (9) in (18) to obtain

$\begin{split}\tilde{H}_{J}(t)=\underbrace{\frac{2(\mu-1)^{\mu/2}(m_{J}(t))^{\mu/2}\,\mathbb{K}_{\mu}(2\sqrt{(\mu-1)m_{J}(t)})}{\Gamma(\mu)}}_{\coloneqq
C_{\mu}(m_{J})},\end{split}$

where $\mathbb{K}_{a}(x)$ stands for the modified Bessel function of the
second kind of order $a$ evaluated at $x$. Hence, if we consider system (10)
under the assumption of weak competence-based contact function we obtain

$\left\\{\begin{aligned}
\frac{\mathrm{d}S(t)}{\mathrm{d}t}&=-\beta\tilde{H}_{S}(t)S(t)\tilde{H}_{I}(t)I(t)+(1-\alpha)\gamma
I(t),\\\
\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\beta\tilde{H}_{S}(t)S(t)\tilde{H}_{I}(t)I(t)-\delta
E(t),\\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t}&=\delta(1-\eta)E(t)-\gamma
I(t),\\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t}&=\delta\eta E(t)+\alpha\gamma
I(t),\end{aligned}\right.$ (19)

which becomes

$\left\\{\begin{aligned} \frac{\mathrm{d}S(t)}{\mathrm{d}t}&=-\beta
C_{\mu}(m_{S})\,C_{\mu}(m_{I})S(t)I(t)+(1-\alpha)\gamma I(t),\\\
\frac{\mathrm{d}E(t)}{\mathrm{d}t}&=\beta
C_{\mu}(m_{S})\,C_{\mu}(m_{I})S(t)I(t)-\delta E(t),\\\
\frac{\mathrm{d}I(t)}{\mathrm{d}t}&=\delta(1-\eta)E(t)-\gamma I(t),\\\
\frac{\mathrm{d}R(t)}{\mathrm{d}t}&=\delta\eta E(t)+\alpha\gamma
I(t).\end{aligned}\right.$ (20)

The next equation will help to close the system

$\int_{\mathbb{R}^{+}}x\,f_{J}^{\infty}(t)\mathrm{e}^{-k/x}\mathrm{e}^{-x}\,\mathrm{d}x=m_{J}(t)J(t)C_{\mu-1}(m_{J}),$
(21)

which is a straightforward consequence of the following property of the
modified Bessel functions of the second kind

$x\,\mathbb{K}_{\mu+1}(x)-2\mu\mathbb{K}_{\mu}(x)=x\,\mathbb{K}_{\mu-1}(x),\quad
x\in\mathbb{R}^{+}.$

Again under the assumptions that $\lambda_{J}=\lambda_{CJ}=\lambda$ for
$J\in\mathcal{C}$, integrating with respect to $x$ equation (11), with the aid
of equation (21), we get

$\left\\{\begin{aligned} \frac{\mathrm{d}[m_{S}(t)S(t)]}{\mathrm{d}t}&=-\beta
C_{\mu-1,S}C_{\mu,k}S(t)m_{S}(t)I(t)+(1-\alpha)\gamma m_{I}(t)I(t)\\\
&\phantom{{}={}}+\lambda(m_{B}-m_{S}(t))S(t),\\\
\frac{\mathrm{d}[m_{E}(t)E(t)]}{\mathrm{d}t}&=\beta
C_{\mu-1,S}C_{\mu,k}S(t)m_{S}(t)I(t)-\delta m_{E}(t)E(t)\\\
&\phantom{{}={}}+\lambda(m_{B}-m_{E}(t))E(t),\\\
\frac{\mathrm{d}[m_{I}(t)I(t)]}{\mathrm{d}t}&=\delta(1-\eta)m_{E}(t)E(t)-\gamma
m_{I}(t)I(t)+\lambda(m_{B}-m_{I}(t))I(t),\\\
\frac{\mathrm{d}[m_{R}(t)R(t)]}{\mathrm{d}t}&=\delta\eta
m_{E}(t)E(t)+\alpha\gamma
m_{I}(t)I(t)+\lambda(m_{B}-m_{R}(t))R(t),\end{aligned}\right.$ (22)

which, using again the fact that

$\frac{\mathrm{d}[m_{J}(t)J(t)]}{\mathrm{d}t}=J(t)\frac{\mathrm{d}m_{J}(t)}{\mathrm{d}t}+m_{J}(t)\frac{\mathrm{d}J(t)}{\mathrm{d}t},$

leads to

$\left\\{\begin{aligned} \frac{\mathrm{d}m_{S}(t)}{\mathrm{d}t}&=-\beta
C_{\mu,I}\,m_{S}(t)I(t)\bigl{[}C_{\mu-1,S}\,m_{S}-C_{\mu,S}\bigr{]}+(1-\alpha)\gamma\frac{I(t)}{S(t)}\bigl{[}m_{I}(t)-m_{S}(t)\bigr{]}\\\
&\phantom{{}={}}+\lambda(m_{B}-m_{S}(t)),\\\
\frac{\mathrm{d}m_{E}(t)}{\mathrm{d}t}&=\beta
C_{\mu,I}\,\frac{I(t)S(t)}{E(t)}\bigl{[}C_{\mu-1,S}\,m_{S}(t)-C_{\mu,S}\,m_{E}(t)\bigr{]}+\lambda(m_{B}-m_{E}(t)),\\\
\frac{\mathrm{d}m_{I}(t)}{\mathrm{d}t}&=\delta(1-\eta)\frac{E(t)}{I(t)}(m_{E}(t)-m_{I}(t))+\lambda(m_{B}-m_{I}(t)),\\\
\frac{\mathrm{d}m_{R}(t)}{\mathrm{d}t}&=\delta\eta\frac{E(t)}{R(t)}(m_{E}(t)-m_{R}(t))+\alpha\gamma\frac{I(t)}{R(t)}(m_{I}(t)-m_{R}(t))+\lambda(m_{B}-m_{R}(t)).\end{aligned}\right.$
(23)

##  4 Examples and applications

In this section we numerically validate the modeling framework proposed in (3)
with local incidence rate 4 in the settings 1-(2). We stress that those form
of contact functions generate different macroscopic models that have been
defined in (15),(17) and (20),(23), respectively, for $\varepsilon\ll 1$. Once
established the consistency of the approach, we proceed by exploiting the
macroscopic sets of equations for calibration purposes based on a freely
available repository for the spreading of hashtags linked to known fake news.
The proposed data-oriented approach is fundamental to experimentally observe
the different impact of the contact function in identifying impact of
competence in the fake news dynamics.

From the numerical point of view we will exploit an implicit structure
preserving method for the Fokker-Planck operator (9) based on the schemes
presented in [32]. The advantage of these methods relies on an arbitrarily
accurate description of the steady state distribution of the Fokker-Planck
model of interest. Similar approaches have been investigated in a different
context also in [12, 13, 30].

###  4.1 Test 1: Validation of the social closure

In this first test we compare the evolution of mass fractions $J(t)$ and means
$m_{J}(t)$, $J\in\mathcal{C}$, obtained from direct integration of
$f_{J}(x,t)$, solution to (3), with respect to the competence
$x\in\mathbb{R}_{+}$, with the macroscopic models (15),(17) and (20)–(23) for
several regimes of $\varepsilon>0$. We start by outlining the procedure by
which we solve the system of kinetic equations (3) with Fokker-Planck
interaction operators. Since $\varepsilon>0$ is assumed to be small, we adopt
a time splitting procedure. In particular, upon introducing a time
discretization $t^{n}=n\Delta t$, $\Delta t>0$ constant, we proceed as
follows.

##### I. Fokker-Planck solver.

At time $t=t^{n}$, we determine the distributions $f_{H}(x,t)$ for all
$H\in\mathcal{C}$ solution to

$\left\\{\begin{aligned}
\partial_{t}F_{H}(x,t)&=\frac{1}{\varepsilon}\mathcal{Q}(F_{H})(x,t),\quad
t\in(t^{n},t^{n+1/2}]\\\
F_{H}(x,t^{n})&=F_{H}^{0}(t^{n},x),\end{aligned}\right.$

where $\mathcal{Q}(f_{H})$ is the Fokker-Planck operator defined in Section
2.3 whose form, in the hypothesis $\lambda_{H}=\lambda_{CJ}=\lambda$, is given
by

$\mathcal{Q}(F_{H})(x,t)=\partial_{x}\left[\lambda(x-m)F_{H}+\dfrac{\sigma^{2}}{2}\partial_{x}(x^{2}F_{H}(x,t))\right].$

In this step we take advantage of an implicit structure preserving (SP) scheme
for Fokker-Planck equations [32] and describes with arbitrary accuracy the
steady state of the model. In Figure 2 we report for several
$\varepsilon=1,10^{-1},10^{-2},10^{-4}$ the numerical solution of the Fokker-
Planck model in the time interval $[0,T]$, $T=2$, obtained from a
discretization of the domain $[0,4]$ with $N_{v}=201$ grid points and with
$\Delta t=10^{-4}$. We may observe that the scheme is capable to approximate
the inverse Gamma analytical equilibrium $f^{\infty}(x)$. We also report the
evolution of the $L^{2}$ numerical error computed as
$\|f(x,t)-f^{\infty}(x)\|_{L^{2}}$ in the time frame $[0,2]$ from which we can
observe how for sufficiently small values of $\varepsilon$ we correctly
approximate the given equilibrium distribution.

Figure 2: Test 1. Left: numerical distribution obtained with SP implicit
scheme at time $T=2$ for several $\varepsilon=1$, $10^{-1}$, $10^{-2}$,
$10^{-4}$. Right: evolution of the $L^{2}$ error
$\|f(x,T)-f^{\infty}(x)\|_{L^{2}}$. In all the tests we considered a
discretization of the domain $[0,4]$ obtained with $N_{v}=201$ grid points and
$\Delta t=10^{-4}$.

##### II. Advection-Reaction step.

Hence, we consider the distribution obtained in the interaction step as an
input for the advection-reaction dynamics for $t\in[t^{n+1/2},t^{n+1}]$

$\left\\{\begin{aligned}
\partial_{t}f_{S}(x,t)&=-K(f_{S},f_{I})(x,t)+(1-\alpha)\gamma
f_{I}(x,t)+\lambda\partial_{x}\left[\left(x-m_{B}\right)f_{S}(x,t)\right],\\\
\partial_{t}f_{E}(x,t)&=K(f_{S},f_{I})(x,t)-\delta
f_{E}(x,t)+\lambda\partial_{x}\left[\left(x-m_{B}\right)f_{E}(x,t)\right],\\\
\partial_{t}f_{I}(x,t)&=\delta(1-\eta)f_{E}(x,t)-\gamma
f_{I}(x,t)+\lambda\partial_{x}\left[\left(x-m_{B}\right)f_{I}(x,t)\right],\\\
\partial_{t}f_{R}(x,t)&=\delta\eta f_{E}(x,t)+\alpha\gamma
f_{I}(x,t)+\lambda\partial_{x}\left[\left(x-m_{B}\right)f_{R}(x,t)\right]\\\
f_{H}(x,t^{n+1/2})&=F_{H}(x,t^{n+1/2}),\quad
H\in\\{S,E,I,R\\}.\end{aligned}\right.$

In particular, we adopted a second order Lax-Wendroff scheme coupled with an
explicit time integration.

In the test of this subsection, unless otherwise specified, we prescribe as
initial datum the distribution

$f_{H}(x,0)=H(t)\dfrac{a_{1}^{a_{2}}}{\Gamma(a_{2})}x^{-a_{2}-1}\exp\\{-a_{1}/x\\},\qquad
H\in\\{S,E,I,R\\},$ (24)

where $a_{1}=2(a_{2}-1)$ and $a_{2}=1.25$ with initial mass fractions

$S(0)=0.98,\qquad E(0)=0.018,\qquad I(0)=R(0)=0.001.$ (25)

We consider the choice of parameters $m=\sum_{H\in\mathcal{C}}H(0)m_{H}(0)$,
$\lambda=0.25$ and $\sigma=0.01$ for (4.1). The fake news dynamics is
regulated by the following choice of parameters $\alpha=0.9$, $\beta=20$,
$\gamma=0.2$, and $\delta=0.05$. For contact rates 1–2 we compared the
evolution of mass fractions and mean values obtained from the integration of
(3) with the ones derived in Section 3. We consider the time interval $[0,T]$,
$T=12$, a uniform time discretization with $\Delta t=10^{-4}$ and
$\varepsilon=1$, $10^{-4}$. In particular, Figure 3 refers to the case
$\kappa(x,x_{*})=\beta/xx_{*}$ and Figure 4 to the case $\kappa(x,x_{*})=\beta
e^{-x-x_{*}}$. In both cases we may observe that for small values of
$\varepsilon$ the obtained macroscopic models are accurate in describing the
trends of observable quantities of the kinetic field model. The macroscopic
systems of coupled ODEs has been solved through a RK4 numerical scheme with
$\Delta t=10^{-4}$.

Figure 3: Test 1. Evolution of mass fractions (left) and mean values (right)
obtained from direct integration of the kinetic model (3) in the case
$\kappa(x,x_{*})=\beta/(x\,x_{*})$, for $\varepsilon=1$, $\varepsilon=10^{-4}$
together with the evolution of mass fractions of the macroscopic model
(15)–(17). In both cases we considered $\alpha=0.9$, $\beta=20$, $\gamma=0.2$,
$\delta=0.05$. The kinetic model has been solved through the scheme
$\mathrm{I}$–$\mathrm{II}$ over the domain $[0,4]$, discretization obtained
with $N_{v}=201$ grid points and $\Delta t=10^{-4}$. The initial distribution
of the kinetic model has been defined in (24)-(25).

Figure 4: Test 1. Evolution of mass fractions obtained from direct integration
of the kinetic model (3) in the case $\kappa(x,x_{*})=\beta e^{-x-x_{*}}$, for
$\varepsilon=1$, $\varepsilon=10^{-4}$ together with the evolution of mass
fractions of the macroscopic model (15)–(17). In both cases we considered
$\alpha=0.9$, $\beta=20$, $\gamma=0.2$, $\delta=0.05$. The kinetic model has
been solved through the scheme $\mathrm{I}$–$\mathrm{II}$ over the domain
$[0,4]$, discretization obtained with $N_{v}=201$ grid points and $\Delta
t=10^{-4}$. The initial distribution of the kinetic model has been defined in
(24)-(25).

###  4.2 Test 2: A data driven application to Twitter

In this test we focus on the spreading of the fake news by considering
available Twitter data from the repository TweetSets111Justin Littman. (2018).
TweetSets. Zenodo. https://doi.org/10.5281/zenodo.1289426. In details, we
analyzed the evolution from March to November, 2020 of the hashtag #facemask
related to the COVID-19 pandemic, and of the hashtags
#hurricaneflorence#fakenews both associated to the hurricane Florence of
September 2018 that caused catastrophic damages in USA, particularly in the
states of North Carolina and South Carolina.

In the following we will assume that the competence variable is strongly
related to the education level of a country. The data for the initial
distribution of education has been extrapolated by the available Italian data
from 2011 ISTAT census, and has been considered as representative data of a
prototypical Western country [21]. As underlined in [21] the cumulative
distribution of education exhibits a power-law type of tail. For this reason,
as an approximation of the competence distribution we considered an inverse
Gamma of the form

$g(x)=\dfrac{c_{1}^{c_{2}}}{\Gamma(c_{2})}\dfrac{e^{-c_{1}/x}}{x^{1+c_{2}}},$
(26)

with $c_{1},c_{2}>0$ obtained by data fitting. More precisely, we measure the
education level on the scale $[0,6]$ where 6 represents the education of
people with a PhD (see Figure 5).

Figure 5: Test 2. Competence distribution and its inverse Gamma approximation
$f(x)$ (26) corresponding to $c_{1}\approx 0.75$, $c_{2}\approx 1.25$ and
leading to a mean competence background of $m_{B}=3$. Data refers to 2011
Italian census and are used as representative of a prototypical Western
country. On the $x$-axis we indicated with $(1)$ lower secondary education,
$(2)$ upper secondary education, $(3)$ undergraduate, $(4)$ master, $(5)$
second level master, $(6)$ doctorate.

####  4.2.1 Test 2A: Fitting the model to data

Once we have obtained the initial competence distribution together with the
value of $m_{B}$ we can estimate the parameters of the models defined in
(15)–(17) and (20)–(23). Several approaches have been proposed in the
literature, see e.g., [25]. It is worth to mention that several uncertainties
are present in data linked to news-monitoring. For example the total
population size is generally unknown and the total number of Twitter accounts
represent an upper bound over the real active users.

The approach adopted in [16], and subsequently in [17, 26], is to treat this
quantity as a parameter to be determined in the minimization process along
with the parameters of the models. To reduce the number of parameters to
optimize we follow a different path. In particular, as initial guess on the
total population size, since the datasets that we used for the fitting were
based on U.S. hashtags, we considered that each fake-news spreader has in
average 453 followers 222Kickfactory, link at
https://kickfactory.com/blog/average-twitter-followers-updated-2016/, last
accessed: 21th February 2022.. Hence, in average we may expect that the total
number of susceptible is given by the total number of tweets multiplied by the
average number of followers. To take also into account both the number of bots
on Twitter as found in [35] (and references therein) and users whose activity
could be not assiduous enough to matter during the lifespan of the considered
fake news, the initial guess was also reduced by a factor of $4$.

Let us denote by $\hat{I}(t)$ the number of active spreaders obtained from the
data, while $I(t)$ is the number of infectious agents given by the macroscopic
differential model. Hence, we consider the following cost functional

$F(\hat{I},I)=\bigg{\lVert}\int_{t_{0}}^{t_{f}}\hat{I}(t)dt-\int_{t_{0}}^{t_{f}}I(t)dt\bigg{\rVert}_{L^{2}([t_{0},t_{f}])},$

where $[t_{0},t_{f}]$ is the time-frame (in hours) during which we solve the
minimization problem

$\min_{\alpha,\beta,\gamma,\delta\in\mathbb{R}_{+}}F(\hat{I},I),$ (27)

whereas $\eta$ was kept fixed and equal to $0.5$.

Since data for the evolution of compartments $S,E,R$ are not at our disposal,
as well is not the initial means value for any of the compartments, we solved
the ODE model on $[t_{\star},t_{0}]$, where $t_{0}$ is the starting point of
the spreading process and $t_{\star}$ is a suitable unknown time previous to
$t_{0}$ starting from single exposed, infectious and recovered individuals.
The idea is to simulate an initial situation for the spread of fake news to
happen. Furthermore, we considered initial mean values equal to the half of
the mean background distribution of competence, i.e. $m_{J}(0)=1.5$.

Figure 6: Test 2A. Top row: optimization results for #facemask. Bottom row: optimization results for #florence#fakenews. From right to left: approximation of raw data on the number of tweets, cumulative distribution and evolution of the mean competence for each compartment.  Parameter | #facemask | #florence#fakenews
---|---|---
| $\kappa(x,x_{*})=\beta/(x\,x_{*})$ | $\kappa(x,x_{*})=\beta e^{-x-x_{*}}$ | $\kappa(x,x_{*})=\beta/(x\,x_{*})$ | $\kappa(x,x_{*})=\beta e^{-x-x_{*}}$
$\alpha$ | 0.9995 | 0.9993 | 1.0000 | 0.9999
$\beta$ | 0.0122 | 0.2937 | 0.0901 | 0.9999
$\delta$ | 0.0237 | 0.0336 | 1.0000 | 0.1930
$\gamma$ | 0.0046 | 0.0079 | 0.2127 | 0.9999
Table 2: Test 2A. Estimated parameters for the entire datasets for the
hashtags #facemask (second and third column) and #florence#fakenews (fourth
and fifth column).

In Figure 6 we compare the evolution on the number of tweets regarding the
hashtag #facemask, from 3rd March 2020 to 22 November 2020, and the hashtag
#florence#fakenews, from 11th September 2018 to 4th October 2018, with the
evolution of the model (15), (17). The obtained parameters are reported in
Table 2.

In both cases, we may observe that the evolution of the mean competence levels
are different in the four compartments and, in particular, that low competence
levels are associated to exposed and infectious agents, i.e., the active
spreaders. The outcome reflects the intuitive idea that the disinformation
could be driven by the lack of capability to recognize an information as
purposely false in the first place.

To better take into account the impact of a competence-based contact rate
function $\kappa(x,x_{*})$, we also computed the associated basic reproduction
number $R_{t}$ using the parameters $\beta$ and $\gamma$ estimated previously
for both datasets, reported in Table 2. Following [4, 19], and omitting the
details for brevity, we consider a generalized version of the classical
reproduction number defined as

$R(t)=\frac{\int_{\mathbb{R}^{+}}K(f_{S},f_{I})(x,t)\,dx}{\int_{\mathbb{R}^{+}}\gamma
f_{I}(x,t)\,dx},$ (28)

where again we leveraged the structure preserving scheme proposed in [32] to
perform the calculations.

Figure 7: Test 2A. Evolution of $R_{t}$ in the first 24 hours of datasets
#facemask (left) and #florence#fakenews (right) for the parameters estimated
in Table 2 relative to the introduced contact functions $\kappa(x,x_{*})$.

####  4.2.2 Test 2B: Forecasting under data uncertainties

To analyze the impact of uncertainties in data and parameters we consider a 3D
random variable $\mathbf{z}=(z_{1},z_{2},z_{3})$ with distribution
$\rho(\mathbf{z})$. We will suppose that the random vector $\mathbf{z}$ has
independent components, i.e.
$\rho(\mathbf{z})=\rho_{1}(z_{1})\rho_{2}(z_{2})\rho_{3}(z_{3})$. Taking into
account parametric uncertainties, we consider the estimated model parameters
as follows

$\beta(z_{1})=\beta_{0}(1+c_{\beta}z_{1}),\qquad\gamma(z_{2})=\gamma_{0}(1+c_{\gamma}z_{2}),\qquad\delta(z_{3})=\beta_{0}(1+c_{\delta}z_{3}),$
(29)

where we supposed $z_{1},z_{2},z_{3}\sim\mathcal{U}([-1,1])$ and
$c_{\beta},c_{\gamma},c_{\delta}>0$. As a result, the macroscopic quantities
describing the evolution of compartments result affected by the introduced
uncertainties increasing their dimensionality $J(\mathbf{z},t)$,
$m_{J}(\mathbf{z},t)$, $J\in\mathcal{C}$. In order to handle efficiently the
introduced uncertainties in the dynamics we adopt a stochastic collocation
approach based on stochastic Galerkin methods, we refer the interested reader
to [38] for an introduction and to [2, 39] for applications in compartmental
modelling of epidemic dynamics. This class of methods allows to accurately
quantify the propagation stochasticity in a parametric differential model when
information on the uncertainties’ distribution are available. We remark that
fast convergence properties hold under suitable regularity assumptions on the
problem’s solution. In details, we construct a 3D sample
$\\{z_{i,k}\\}_{k=0}^{M}$, $i=1,2,3$, obtained in a collocation setting
through Gauss-Legendre polynomials with $M=5$ nodes.

In Figure 8 we display the dynamics of the considered fake-news with respect
to available data. In details, for #florence#fakenews we consider the period
from September 11th 2018 to September 21st 2018. We consider two successive
prediction horizons respectively of 1 day, i.e. the parameters of the models
are calibrated taking into account data until September 19th, and a 2 days
prediction horizon, where the calibration is based only on data until
September 18th. Regarding #facemask we considered the period from March 3rd to
May 17th. Also in this case we consider two successive prediction horizons of
1 week, i.e. the parameters of the models are calibrated taking into account
data until May3rd, and a two weeks prediction horizon, where the calibration
is based on data until May10th.

We highlight in dashed black and magenta the expected value of the predicted
number of tweets
$\mathbb{E}[I(\mathbf{z},t)]=\int_{0}^{1}I(\mathbf{z},t)\rho_{1}(z_{1})\rho_{2}(z_{2})\rho_{3}(z_{3})dz_{1}dz_{2}dz_{3}$.
Together with the expected trends we plot the $95\%$ confidence intervals (CI)
with respect to the random parameters $\beta(z_{1})$, $\gamma(z_{2})$, and
$\delta(z_{3})$. The blue shaded band is relative to the variability in
$\gamma(z_{2})$, the green shaded to the variability in $\delta(z_{3})$
whereas the shaded red is relative to the variability in $\beta(z_{1})$.

Figure 8: Test 2B. Top row: comparison between 14 days (May 3rd–May 14th,
2020) and 7 days (May 3rd–May 10th, 2020) predictions of #facemask based on
the model (15), (17) with uncertain parameter (29) with
$c_{\beta}=200\beta_{0}$, $c_{\gamma}=2\gamma_{0}$, $c_{\delta}=2\times
10^{4}\delta_{0}$. Bottom row: comparison between 24 hours (September 18th-
September 19th, 2018) and 48 hours (September 18th-September 20th, 2018)
predictions of #florence#fakenews based on the model (15), (17) with uncertain
parameter (29) with $c_{\beta}=0.2\beta_{0}$, $c_{\gamma}=\gamma_{0}$,
$c_{\delta}=2\times 10^{4}\delta_{0}$.

###  4.3 Test 3: Competence background in misinformation

In this test we perform a retrospective analysis to study how the background
could influence the dissemination of fake news as a result of a different
learning process. We recall that the background modifies through a learning
dynamic the effectiveness of the level of knowledge in identifying fake news.
As a consequence high values of the background correspond to a high level of
effectiveness of the competence while low values will make it difficult to
identify the fake-news. Indirectly, the background acts as a control term
which limits the spread of the misinformation. This can also be interpreted as
a process of education specific to the identification of fake news that allows
to limit the so-called knowledge neglect phenomenon [18].

We consider the two datasets for the hashtags #facemask and #florence#fakenews
with the estimated parameters reported in Table 2 and we increase the value of
the competence level attained by the background, i.e., $m_{B}$, while keeping
fixed the parameters during the dynamics defined by (15),(17) and (20),(23).

Hence, we performed the test with both choices of a strong and weak competence
based contact function; the results are summarized in Figure 9. In all cases,
we see how increasing the competence of the background reduces the spread of
fake news, leading to a decrease in the cumulative number of tweets of
infectious agents proportional to the increase in the value of $m_{B}$.
Indeed, we can observe how increasing the competence of the background, we
obtain an evident decrease in the overall misinformation for both the examples
considered #facemask and #florence#fakenews.

Figure 9: Test 3. Total number of infectious agents for the hashtags #facemask
(left) and #florence#fakenews (right) as a function of the competence
background. In both cases, we employed the parameters reported in Table 2.

## Concluding remarks

Despite the digital transformation of governments and the modernization of
public administration, a global decline in democracy is occurring around the
world. The spread of fake news created for the purpose of polarizing society
in certain directions poses a risk to democratic institutions. The role of
individuals’ knowledge and the ability to use it in identifying false
information is deemed of paramount importance.

In this paper starting from a model for the description of fake-news
dissemination in the presence of heterogeneous agents with different levels of
competence, through the tools of kinetic theory, reduced-order models have
been derived that allow to keep the effects of the of competence in the
dynamics and that, thanks to their simplified structure, can be interfaced
with data.

The starting model is inspired, as in much of the literature related to fake-
news, to the epidemiology, so it is based on a compartmental structure. The
introduction of competence allows to analyze complex phenomena of great
relevance in contemporary society, such as the effectiveness of control
actions taken to limit the spread of fake news and the role of knowledge
neglect in misinformation.

The methodology adopted in this article is fully general and depends closely
on the equilibrium state of the social variable and the social interaction
function at the basis of fake-news spreading. As a consequence, additional
social variables that play a key role in the spread of misinformation may be
embedded in the dynamics using similar arguments. The ability to have a model
that can be interfaced with the available data allowed us to present some
preliminary examples of applications to the case of fake-news spreading on
Twitter.

## Competing interests

The authors declare no competing interests.

## Data availability statement

The datasets generated during the current study is available from the
corresponding author on reasonable request. The datasets analysed during the
current study are freely available from the website
https://doi.org/10.5281/zenodo.1289426.

### Acknowledgments

This work has been written within the activities of the GNFM and GNCS groups
of INdAM (National Institute of High Mathematics). M.Z. acknowledge partial
support of MUR-PRIN2020 Project “Integrated mathematical approaches to socio-
epidemiological dynamics”. The research of M.Z. was partially supported by
MIUR, Dipartimenti di Eccellenza Program (2018–2022), and Department of
Mathematics “F. Casorati”, University of Pavia. The research of L.P. was
partially supported by FIR project “No hesitation. For effective communication
of Covid-19 vaccination”, University of Ferrara.

## References

  * [1] G. Albi, G. Bertaglia, W. Boscheri, G. Dimarco, L. Pareschi, G. Toscani, M. Zanella.  Kinetic modelling of epidemic dynamics: social contacts, control with uncertain data, and multiscale spatial dynamics.  To appear in: _Predicting Pandemics in a Globally Connected World_ , Vol. 1, N. Bellomo and M. Chaplain Editors, Springer-Nature (2022).
  * [2] G. Albi, L. Pareschi, M. Zanella. Control with uncertain data of socially structured compartmental epidemic models. J Math Biol. 82(7):63.
  * [3] G. Albi, L. Pareschi, M. Zanella. Opinion dynamics over complex networks: kinetic modelling and numerical methods. _Kinet. Relat. Models_ , 10(1): 1-32, 2017
  * [4] G. Bertaglia and L. Pareschi. Hyperbolic compartmental models for epidemic spread on networks with uncertain data: Application to the emergence of Covid-19 in Italy. Math. Mod. Meth. Appl. Sci., 31(12):2495–2531, 2021.
  * [5] L. Boudin, L. Trussardi, Concentration effects in a kinetic model with wealth and knowledge exchanges, Preprint _hal-02470191_ , 2021.
  * [6] C. Cercignani. _The Boltzmann Equation and its Applications_ , Springer: Berlin, Germany, 1988.
  * [7] S.H. Choi, H. Seo. Rumor spreading dynamics with an online reservoir and its asymptotic stability. Netw. Heterog. Media 16(4): 535–552, 2021.
  * [8] S. Cordier, L. Pareschi, and G. Toscani.  On a kinetic model for a simple market economy.  J. Stat. Phys., 120(1):253–277, 2005.
  * [9] E. Cristiani, A. Tosin. Reducing complexity of multiagent systems with symmetry breaking: an application to opinion dynamics with polls. _Multiscale Model. Simul._ , 16(1):528-549, 2018.
  * [10] A. Cossard, G. De Francisci Morales, K. Kalimeri, Y. Mejova, D- Paolotti, M. Starnini. Falling into the Echo Chamber: The Italian Vaccination Debate on Twitter. Proceedings of the International AAAI Conference on Web and Social Media, 14(1):130-140, 2020.
  * [11] D. J. Daley, D. G. Kendall. Epidemics and rumors. _Nature_ , 204:1118, 1964.
  * [12] G. Dimarco, B. Perthame, G. Toscani, M. Zanella. Kinetic models for epidemic dynamics with social heterogeneity. J. Math. Biol., 83:4, 2021.
  * [13] G. Dimarco, G. Toscani, M. Zanella. Optimal control of epidemic spreading in presence of social heterogeneity. _Phil. Trans. R. Soc. A_ , 380:20210160, 2022.
  * [14] B. Düring, M.-T. Wolfram. Opinion dynamics: inhomogeneous Boltzmann-type equations modelling opinion leadership and political segregation. _Proc. R. Soc. A_ , 471: 20150345, 2015.
  * [15] B. Düring, O. Wright. On a kinetic opinion formation model for pre-election polling. _Phil. Trans. R. Soc. A_ , 380:20210154, 2022.
  * [16] J. Fang, E. Dougherty, P. Saraf, C. Yang, N. Ramakrishnan. Epidemiological Modeling of News and Rumors on Twitter. _SNAKDD ’13: Proceedings of the 7 th Workshop on Social Network Mining and Analysis_, 8: 1–9, 2013.
  * [17] J. Fang, W. Wei, Z. Liang, E. Dougherty, C. Yang, L. Chang-Thien, N. Ramakrishnan. Misinformation Propagation in the Age of Twitter. Computer, 47(12):90–94, 2014.
  * [18] L.K. Fazio, N.M. Brashier, B.K. Payne, E. Marsh. Knowledge does not protect against illusory truth. _J. Journal of Experimental Psychology: General_ , 144(5):993–1002, 2015.
  * [19] J. Franceschi, L. Pareschi. Spreading of fake news, competence, and learning: kinetic modeling and numerical approximation, _Phil. Trans. R. Soc. A_ , 380:20210159, 2022.
  * [20] P. Glister. _Digital Literacy_. Wiley Computer Pub., New York, 1997.
  * [21] S. Gualandi, G. Toscani. Pareto tails in socio-economic phenomena: a kinetic description. Economics: The Open-Access, Open-Assessment E-Journal, 12, 2018.
  * [22] A. M. Guess, M. Lerner, B. Lyons, J. M. Montgomery, B. Nyhan, J. Reifler, N. Sircar. A digital media literacy intervention increases discernment between mainstream and false news in the United States and India. _PNAS_ 117(27):15536–15545, 2020.
  * [23] B.I. Hong, N. Hahm, S.H. Choi. SIR Rumor spreading model with trust rate distribution, Netw. Heterog. Media, 3(3): 515–530, 2018.
  * [24] J.B. Horrigan, Digital readiness gaps. Pew Research Center, October 2019. www.pewresearch.org/internet/2016/09/20/digital-readiness-gaps/
  * [25] Y. Liu, A.A. Gayle, A. Wilder-Smith, J. Rocklöv. The reproductive number of COVID-19 is higher compared to SARS coronavirus. J. Travel Med., 27(2): 1–4, 2020.
  * [26] M. Maleki, E. Mead, M. Arani, N. Agarwal. Using an Epidemiological Model to Study the Spread of Misinformation during the Black Lives Matter Movement. Preprint arXiv:2103.12191, 2021.
  * [27] T. Murayama, S. Wakamiya, E. Aramaki, R. Kobayashi. Modeling the spread of fake news on Twitter. PLoS One, 16(4): 1–16, 2021.
  * [28] L. Pareschi, G. Toscani.  Interacting Multiagent Systems: Kinetic equations and Monte Carlo methods.  Oxford University Press, 2013.
  * [29] L. Pareschi, G. Toscani. Wealth distribution and collective knowledge. A Boltzmann approach. _Phil. Trans. R. Soc. A_ 372: 20130396, 2014.
  * [30] L. Pareschi, T. Trimborn, and M. Zanella. Mean-field control variate methods for kinetic equations with uncertainties and applications to socio-economic sciences. Int. J. Uncertain. Quantif., 12(1):61–84, 2022.
  * [31] L. Pareschi, P. Vellucci, M. Zanella. Kinetic models of collective decision-making in the presence of equality bias. _Phys. A_ 467: 201-217, 2017.
  * [32] L. Pareschi, M. Zanella. Structure preserving schemes for nonlinear Fokker-Planck equations and applications. J. Sci. Comput., 74:1575–1600, 2018.
  * [33] G. Toscani. Kinetic models of opinion formation. Commun. Math. Sci., 4(3):481–496, 2006.
  * [34] G. Toscani, A. Tosin, M. Zanella. Opinion modeling on social media and marketing aspects. _Phys. Rev. E_ , 98(2): 022315, 2018.
  * [35] J. Uyheng, K.M. Carley. Bots and online hate during the COVID-19 pandemic: case studies in the United States and the Philippines. J. Comput. Soc. Sc. 3:445–468, 2020.
  * [36] C. Villani. On a new class of weak solutions to the spatially homogeneous Boltzmann and Landau equations. _Arch. Ration. Mech. Anal._ 143, 273–307,1998.
  * [37] S. Vosoughi, D. Roy, S. Aral. The spread of true and false news online, _Science_ 359(6380):1146–1151, 2018.
  * [38] D. Xiu. _Numerical Methods for Stochastic Computations: A Spectral Method Approach_ , Princeton University Press, 2010.
  * [39] M. Zanella, C. Bardelli, G. Dimarco, S. Deandrea, P. Perotti, M. Azzi, S. Figini, G. Toscani. A data-driven epidemic model with social structure for understanding the COVID-19 infection on a heavily affected Italian Province. _Math. Mod. Meth. Appl. Sci._ , 31(12):2533-2570, 2021.