# The threshold model with anticonformity under random sequential updating

Bartlomiej Nowak, Michel Grabisch, Katarzyna Sznajd-Weron

► To cite this version:

Bartlomiej Nowak, Michel Grabisch, Katarzyna Sznajd-Weron. The threshold model with anticonformity under random sequential updating. 2022. halshs-03573376

HAL Id: halshs-03573376

https://shs.hal.science/halshs-03573376

Submitted on 14 Feb 2022

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
CES

Centre d'Économie de la Sorbonne
UMR 8174

The threshold model with anticonformity
under random sequential updating

Bartlomiej NOWAK, Michel GRABISCH,
Katarzyna SZNAJD-WERON

2022.04
# The threshold model with anticonformity under random sequential updating

Bartłomiej Nowak,¹, * Michel Grabisch,², † and Katarzyna Sznajd-Weron¹, ‡

¹Department of Theoretical Physics, Faculty of Fundamental Problems of Technology, Wrocław University of Science and Technology, 50-370 Wrocław, Poland

²University of Paris I Panthéon-Sorbonne, Paris School of Economics 106-112, Bd de l'Hôpital, Paris, France
(Dated: January 28, 2022)

We study an asymmetric version of the threshold model with anticonformity under asynchronous update mode that mimics continuous time. We study this model on a complete graph using three different approaches: mean-field approximation, Monte Carlo simulation, and the Markov chain approach. The latter approach yields analytical results for arbitrarily small systems, in contrast to the mean-field approach, which is strictly correct only for an infinite system. We show that for sufficiently large systems, all three approaches produce the same results, as expected. We consider two cases: (1) homogeneous, in which all agents have the same tolerance threshold, and (2) heterogeneous, in which the thresholds are given by a beta distribution parametrized by two positive shape parameters $\alpha$ and $\beta$. The heterogeneous case can be treated as a generalized model that reduces to a homogeneous model in special cases. We show that particularly interesting behaviors, including social hysteresis and critical mass, arise only for values of $\alpha$ and $\beta$ that yield the shape of the distribution observed in real social systems.

## I. INTRODUCTION

Within the broad class of two-state dynamics [1], threshold models are particularly useful for describing various social and economical phenomena [2-4]. As other binary-state opinion dynamics [5], the threshold model describes the social influence in decision-making for the choice between precisely two alternatives, often denoted by 1 (agree, adopt, be active, etc.) and 0 (disagree, refuse, be inactive, etc.). Although a binary decision framework seems to be oversimplified, it is relevant to surprisingly many complex problems [3].

In the original threshold models of collective behavior, proposed by Schelling [6] and Granovetter [2], an agent takes action 1 if the proportion of his neighbors in state 1 exceeds some threshold, otherwise action 0 is taken. It means that an agent at state 1 may return to state 0, because not enough neighbors are active. On the other hand, in many other threshold models, the transition from state 1 to state 0 is forbidden [3, 7, 8].

Here we will use the original formulation, in which a transition from 1 to 0 is possible, as in [2, 9, 10], but additionally in the presence of anticonformity. Such a model has been already studied from a mathematical point of view under the synchronous update mode [10]. The study was focused on finding absorbing classes, cycles, etc. In this paper, we investigate the same model but under random sequential updating, which mimics continuous time. Contrarily to [10], we focus on phase transitions and phase diagrams, which is a typical approach for statistical physics of opinion formation [11-15].

We study the model on the complete graph, which enables to obtain exact results within the mean-field ap

proach. Independently, we conduct Monte Carlo simulations to validate the theoretical approach. Finally, we present a Markov Chain approach, which not only allows us to obtain results for arbitrary small systems, but also to derive the stationary distribution of visited states.

## II. MODEL

We consider a society of $n$ agents placed at the vertices of an arbitrary graph $G = (N, E)$, where $N = \{1, ..., n\}$ is a set of vertices (agents) and $E$ is the set of undirected edges. Each agent $i$ has a set of neighbors $K_i = \{j \in N : \{i, j\} \in E\}$, and cardinality of this set $|K_i| = k_i$ is the degree of agent $i$. As in many other models, an agent can be in one of two alternative states: 1 (agree, adopt, be active, etc.) or 0 (disagree, refuse, be inactive, etc.). Following [10], we use the term “active” for agents in state 1, and “inactive” for agents in state 0, and denote by $a_i(t)$ the state (action) taken by agent $i$ at time $t$.

We consider two types of social response, anticonformity and conformity, occurring with complementary probabilities $p$ and $1-p$ respectively. In both cases, an agent can change its state, if the ratio of active neighbors is above its tolerance threshold $r_i \in [0, 1]$. Threshold $r_i$ of each agent is the realization of the random variable $R$ with arbitrary distribution function $F_R(r)$ and does not change in time. In case of conformity, an agent follows the others, whereas in case of anticonformity he takes an opposite state to others. Therefore, the dynamics of the agent's state in case of conformity can be written as [10]:

$$a_i(t + \Delta t) = \begin{cases} 1, & \text{if } \frac{1}{k_i} \sum_{j \in K_i} a_j(t) \ge r_i \\ 0, & \text{otherwise,} \end{cases} \quad (1)$$

<footnote>* bartłomiej.nowak@pwr.edu.pl</footnote>

<footnote>† michel.grabisch@univ-paris1.fr</footnote>

<footnote>‡ katarzyna.weron@pwr.edu.pl</footnote>
whereas in case of anticonformity [10]:

$$a_i(t + \Delta t) = \begin{cases} 0, & \text{if } \frac{1}{k_i} \sum_{j \in K_i} a_j(t) \ge r_i \\ 1, & \text{otherwise} \end{cases} \quad (2)$$

In this paper, we use the random sequential update mode, which means that an elementary update consists of:

(i) random drawing of agent i from all n agents

(ii) with probability p agent i anticonforms to the neighborhood, i.e., takes action $a_i(t+\Delta t)$ according to Eq. (2)

(iii) with complementary probability $1-p$ agent i conforms to the neighborhood, i.e., takes action $a_i(t+\Delta t)$ according to Eq. (1)

(iv) time is updated: $t := t + \Delta t$

As usually, $\Delta t = 1/n$ which means that the time unit consists of n elementary updates, which corresponds to one Monte Carlo step (MCS).

### III. TRANSITION PROBABILITIES

Since we limit our study to the complete graph, we can fully describe the state of the system by a single random variable:

$$c = \frac{n_1}{n}, \quad (3)$$

where $n_1$ is the number of agents in state 1 and thus c is the ratio of active agents. Therefore, there are $n+1$ possible states of the system: $0, \frac{1}{n}, \frac{2}{n}, \dots, 1$.

Because we use the sequential (asynchronous) update mode, at most one agent can change its state at a time, and thus we can introduce the following transition probabilities:

$$\begin{aligned} \gamma^{+}(c) &= \mathrm{Pr} \left( c(t + \Delta t) = c(t) + \frac{1}{n} \right), \\ \gamma^{-}(c) &= \mathrm{Pr} \left( c(t + \Delta t) = c(t) - \frac{1}{n} \right). \end{aligned} \quad (4)$$

For our model, the explicit form of these probabilities can be written, according to the algorithm described in the previous section, as follows:

$$\begin{aligned} \gamma^{+}(c) &= (1-p)(1-c)\mathrm{Pr}(R \le c) + p(1-c)\mathrm{Pr}(R > c), \\ \gamma^{-}(c) &= (1-p)c\mathrm{Pr}(R > c) + pc\mathrm{Pr}(R \le c), \end{aligned} \quad (5)$$

where $Pr(R \le c)$ is the probability that the concentration c of active agents is bigger than or equal to the threshold R of the considered agent. This probability is simply the value of the cumulative distribution function $F_R(r)$ at $r = c$. Similarly, $Pr(R > c)$ is the probability

that the concentration of active voters does not exceed the threshold of considered agents and thus it is equal to $1 - F_R(c)$. Therefore, we obtain

$$\begin{aligned} \gamma^{+}(c) &= (1-p)(1-c)F_{R}(c) + p(1-c)(1-F_{R}(c)), \\ \gamma^{-}(c) &= (1-p)c(1-F_{R}(c)) + pcF_{R}(c). \end{aligned} \quad (6)$$

As can be seen from Eq. (4), the concentration of active agents c is a random variable. However, we can easily write the evolution equation for the expected value of c. Moreover, for $n \to \infty$ we can assume that c localizes to the expectation value. Therefore, we can write [5]:

$$\frac{dc}{dt} = \gamma^{+}(c) - \gamma^{-}(c) \quad (7)$$

As usually, we focus on the steady states, i.e., those for which

$$\frac{dc}{dt} = 0. \quad (8)$$

In the next two sections, we will use the condition (8) to calculate the stationary concentration of active agents for two cases: (1) the homogeneous one, in which all agents have the same tolerance threshold (2) the heterogeneous one, in which the distribution of thresholds $F_R(r)$ is given by the beta distribution. We will compare the analytical results with the results of Monte Carlo simulations for the system of size $n = 10^4$, averaged over 10 independent runs collected after $10^4$ Monte Carlo steps. For the Monte Carlo simulations, two types of initial conditions will be used to reproduce all stable solutions of Eq. (8): (1) all agents initially active, which will be denoted by $c(0) = 1$ and (2) all agents initially inactive, which will be denoted by $c(0) = 0$.

### IV. ONE THRESHOLD

In this case, the random variable R takes one value for all agents in the system, i.e., all voters have the same threshold r:

$$\begin{aligned} F_{R}(c) &= \mathbf{1}_{\{r \le c\}}, \\ 1 - F_{R}(c) &= \mathbf{1}_{\{r > c\}}, \end{aligned} \quad (9)$$

where $\mathbf{1}_{\{r \le c\}} = 1$ when $r \le c$ and 0 otherwise. Inserting (9) to Eq. (6) and then to Eq. (7) we obtain:

$$\begin{aligned} \frac{dc}{dt} = (1-p) &[(1-c)\mathbf{1}_{\{r \le c\}} - c\mathbf{1}_{\{r > c\}}] \\ &+ p[(1-c)\mathbf{1}_{\{r > c\}} - c\mathbf{1}_{\{r \le c\}}] \end{aligned} \quad (10)$$

From (10), we obtain several trivial fixed points:

$$(p = 0, c = 0) \qquad \forall r \ne 0 \quad (11)$$

$$(p = 1, c = 0) \qquad r = 0 \quad (12)$$

$$(p = 0, c = 1) \qquad \forall r \in [0, 1] \quad (13)$$

$$(p = 0.5, c = 0.5) \qquad \forall r \in [0, 1] \quad (14)$$

$$(p = 1 - r, c = r) \qquad \forall r \in [0, 1] \quad (15)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Dependency between the stationary concentration of active agents c and the probability of anticonformity p for model with one threshold for different values of the parameter r (indicated in the plots). Solid lines represent stable fixed points obtained analytically from Eq. (16). Symbols represent Monte Carlo simulations from two initial conditions indicated in the legend.</figcaption></figure>

The remaining solutions can be obtained by solving Eq. (8), which leads to

$$p = \frac{\mathbf{1}_{\{r \le c\}} - c}{\mathbf{1}_{\{r \le c\}} - \mathbf{1}_{\{r > c\}}}, \quad (16)$$

what is equivalent to the following cases

$$\forall c < r \quad p = c, \quad \forall c \ge r \quad p = 1 - c. \quad (17)$$

From the above analysis we do not obtain steady state for any value of $r \in [0, 1]$ if $p > r$. However, from the evolution of Eq. (10), as well as from the Monte Carlo simulations, it seems that the system approaches the state $c = r$ for $p > r$, which is shown in Fig. 1. This raises the question of what the evolution of the system actually looks like and how it can be understood on the basis of transition probabilities $\gamma^+(c), \gamma^-(c)$. The above question can be answered if we split the transition probabilities into cases

$$\begin{aligned} \forall c < r \quad \gamma^+(c) &= p(1-c) \wedge \gamma^-(c) = (1-p)c, \\ \forall c \ge r \quad \gamma^+(c) &= (1-p)(1-c) \wedge \gamma^-(c) = pc. \end{aligned} \quad (18)$$

Doing this we easily observe that they do not cross at any point, when $\forall r > 0.5 \ p \ge r$ or $\forall r \le 0.5 \ p > r$, see the fourth column of Fig. 2. This implies no steady state. On the other hand, for $p < r$ transition probabilities $\gamma^-(c)$ and $\gamma^+(c)$ cross each other, as shown in the first three

columns of Fig. 2, i.e., the steady state $\gamma^-(c) = \gamma^+(c)$ exists.

There is another illustrative way to visualize the dynamics of the system based on the idea of potential $V(c)$ [16]:

$$V(c) = - \int F(c)dc = - \int \frac{dc}{dt}dc, \quad (19)$$

where $F(c)$ plays the role of a generalized force, which drives the dynamics of the system. Using such an approach, we draw a ball sliding down the walls of a potential well [16], as shown in Fig. 2. To calculate the explicit form of $V(c)$ we use Eq. (18), which leads to:

$$\begin{aligned} \forall c < r \quad V(c) &= \frac{c^2}{2} - cp, \\ \forall c \ge r \quad V(c) &= \frac{c^2}{2} - c(1-p). \end{aligned} \quad (20)$$

The steady states are local extrema of $V(c)$. From Eq. (20) we see that the potential has a discontinuity at $c = r$ which implies no maximum (unstable steady state). Still, at most two minima (stable steady states) are possible. In general, the number of minima, denoted by $M(r, p)$, can be described as follows

$$\forall r > 0.5 \quad M(r, p) = \begin{cases} 2, & \text{for } p \le 1 - r \\ 1, & \text{for } 1 - r < p < r \\ 0, & \text{for } p \ge r, \end{cases} \quad (21)$$

$$\forall r \le 0.5 \quad M(r, p) = \begin{cases} 2, & \text{for } p < r \\ 1, & \text{for } r \le p \le 1 - r \\ 0, & \text{for } p > r. \end{cases} \quad (22)$$

In conclusion, despite the lack of steady state in the case $M(r, p) = 0$ we can observe the flow of the system is towards the point $c = r$. It reaches an asymptotic minimum at this point because from both, the left and right boundaries, the system flow is towards this minimum. This explains the behavior shown in Fig. 1, which was at first incomprehensible and inspired the above analysis.

## V. BETA DISTRIBUTION

In the previous section, we studied the homogeneous system, in which all agents had the same value of the tolerance threshold $r$. However, we can also consider more general distributions of thresholds, allowing for heterogeneity. The most useful are distributions whose support values $r \in [0, 1]$ and show a variety of shapes. This is the case of the beta distribution with two parameters $\alpha$ and $\beta$, considered previously, for the models of tolerance without anticonformity [4]. It has a well-defined cumulative distribution function

$$F_R(r) = I_r(\alpha, \beta) = \frac{B(r, \alpha, \beta)}{B(\alpha, \beta)}, \quad (23)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2. Analysis of the steady states and the stability of the system for two values of threshold $r = 0.6$ (two first rows) and $r = 0.5$ (two last rows) and four values of $p = 0.2$ (first column), $p = 0.4$ (second column), $p = 0.45$ (third column), $p = 0.8$ (fourth column). In first and third row, solid lines represents values of $\gamma^+$ and dotted lines stands for $\gamma^-$ obtained with Eqs. (18). Potentials $V(c)$ (second and fourth row) are obtained with Eqs. (20). In all subplots, filled circles denote continuity, while empty circles denote lack of continuity at this point.</figcaption></figure>

where $I_r(\alpha, \beta)$ is the regularized incomplete beta function, which can be defined in terms of the incomplete beta function $B(r, \alpha, \beta)$ and the complete beta function $B(\alpha, \beta)$. Inserting $F_R(r)$ given by Eq. (23) to (6) we obtain the transition probabilities $\gamma^+(c), \gamma^-(c)$. Then inserting them to Eq. (7) we get

$$ \frac{dc}{dt} = (1-p)[(1-c)I_c(\alpha, \beta) - c(1-I_c(\alpha, \beta))] + p[(1-c)(1-I_c(\alpha, \beta)) - cI_c(\alpha, \beta)]. \quad (24) $$

Again, we can point out the obvious steady states ($p = 0, c = 0$), ($p = 0, c = 1$) for arbitrary values of $\alpha$ and $\beta$.

For $c = 1/2$ formula (24) boils down to

$$ \left. \frac{dc}{dt} \right|_{c=\frac{1}{2}} = \left( I_{\frac{1}{2}}(\alpha, \beta) - \frac{1}{2} \right) (1 - 2p), \quad (25) $$

which has two roots. The first one $p = 1/2$ gives the fixed point ($p = \frac{1}{2}, c = \frac{1}{2}$). The other root $I_{\frac{1}{2}}(\alpha, \beta) = \frac{1}{2}$ exists if the beta distribution is symmetric around the value $\frac{1}{2}$. This happens for $\alpha = \beta$, which leads to the conclusion that the value $c = \frac{1}{2}$ is a fixed point for all values of $p$ if $\alpha = \beta$. For all remaining solutions we have the following
relation:

$$p = \frac{I_c(\alpha, \beta) - c}{2I_c(\alpha, \beta) - 1}. \qquad (26)$$

The information about the stability of the steady state is given by the sign of the derivative

$$\frac{dF}{dc} = \frac{c^{\alpha-1}(1-c)^{\beta-1}\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}(1-2p) - 1. \qquad (27)$$

The state is stable if the above derivative is negative and unstable otherwise. The overall behavior of the model is summarized in Fig. 3. In the insets of this figure the dependence between the stationary value of c and parameter p is shown. Two shaded areas in Fig. 3 correspond to the situation in which at least one of the parameters α, β is smaller than 1. In this case, for $p > 0$ there is always only one steady state and c is monotonically increasing ($\beta < \alpha$), monotonically decreasing ($\beta > \alpha$) or constant ($\beta = \alpha$) function of p.

Recalling the shape of the probability density function (PDF) of the beta distribution, we can draw some conclusions. If the PDF of the tolerance threshold is a monotonically decreasing function of the threshold r, then the concentration of active agents decreases with the probability of anticonformity p, and vice versa. If the PDF has the highest values at r = 0 and r = 1, being a convex function of r, then for all values of p > 0 the stationary value of active agents is 0.5.

The most complex behavior is seen if both shape parameters α, β are greater than 1 but not infinitely large, which corresponds to a uni-modal PDF, with zero probabilities at both end of the interval range, i.e., at r = 0 and r = 1. This case correspond to moderate tolerance [4], and it is a typical shape of the distribution of actual trait manifestations in behavior, as reported by psychologists [17]. In such a case, the phase transitions appear, as shown in Fig. 3. As long as $\beta = \alpha$, which corresponds to the symmetric PDF, there is a continuous phase transition between the phase in which one type of agent (active or inactive) dominates, and the symmetrical phase without the domination. The critical point, at which this transition occurs, can be calculated by solving the equation

$$\left. \frac{dF}{dc} \right|_{c=\frac{1}{2}} = 0, \qquad (28)$$

which gives:

$$p_1^* = \frac{1}{2} - \frac{\Gamma^2(\alpha)}{2^{1-2(\alpha-1)}\Gamma(2\alpha)}. \qquad (29)$$

For $\alpha \neq \beta$, as long as shape parameters are finite and at least one of them is larger than 1, we obtain an interesting behavior, with the jump at some value of $p = p^*$ and hysteresis, as shown in Fig. 3. This can be especially useful to describe the innovation diffusion. For example,

if $\beta > \alpha$ then for the small value of $p < p^*$ there is possibility of high adoption if the initial fraction of adopted is above the critical mass. However, if the initial fraction of adopted is too low, i.e. below the critical mass, the innovation cannot spread in the society. Similar behavior has been recently reported for the completely different mathematical model of the collective decision-making with social learners for unequal merit options [18].

It is worth noticing that for $\alpha, \beta \to \infty$ we can recover the solution for the model with one threshold, as shown in Fig. 3. We are able to do that by recalling the formula giving the mode of the beta distribution with $\alpha, \beta > 1$:

$$m = \frac{\alpha - 1}{\alpha + \beta - 2}. \qquad (30)$$

While $\alpha, \beta \to \infty$, the beta distribution is a 1-point degenerate distribution with probability 1 at the midpoint m and 0 elsewhere. Thus, to obtain the case with the mode at the point m = r, i.e., recover the distribution for one threshold, parameters α and β should follow the formula

$$\beta = \frac{(1 - r)\alpha - 1 + 2r}{r} \qquad (31)$$

for $\alpha, \beta \to \infty$.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Phase diagram for the heterogeneous model with thresholds described by the beta distribution parametrized by two shape parameters α and β. Each inset shows representative behavior of c(p) for a given area of the phase diagram. Solid lines in the insets correspond to stable stationary states, whereas dashed lines correspond to unstable stationary states.</figcaption></figure>

All results obtained analytically for beta distribution can be also obtained by Monte Carlo simulations, as shown in Fig. 4.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Representative dependencies between the stationary concentration of spins up and the probability of anticonformity for model with beta distribution for different values of the parameters α and β: (a) α = β ≤ 1, (b) α = β > 1, (c) α < β < 1; α close to β, (d) 1 > α > β; α close to β, (e) α ≤ 1 ∧ α < β, (f) β ≤ 1 ∧ α > β, (g) 1 < α < β, (h) α > β > 1. Solid and dotted lines represent stable and unstable steady states respectively, obtained with Eq. (26). The exact values of parameters in the plots are as follows: (a) α = β = 0.9, (b) α = β = 4, (c) α = 0.1; β = 0.2, (d) α = 0.2; β = 0.1, (e) α = 1; β = 3, (f) α = 3; β = 1, (g) α = 5; β = 8, (h) α = 8; β = 5. Symbols represent Monte Carlo simulations from two initial conditions, denoted in the legend. The results are averaged over 10 runs and collected after 10⁴ MCS for system of size 10⁴.</figcaption></figure>

VI. MARKOV CHAIN APPROACH

Previously, we were assuming that the size of the sys-
tem is infinite, i.e., $n \to \infty$. However, such an assumption
is not very realistic for social systems. Actually, social
scientists are often interested in small systems. There-
fore, in this section, we make analysis of the convergence
of $c$ in the long run using Markov chains for arbitrary

small systems. The advantage of the Markov chain ap-
proach in the context of agent-based modeling of opinion
dynamics has been already reported in [19].

Transition probabilities given by Eq. (4) allows us to
write the $(n+1) \times (n+1)$ transition matrix, whose general
term $(i,j)$ indicates the probability of transition from
state $i$ to state $j$. Due to the asynchronous update mode,
$\mathbf{P}$ is a tridiagonal row-stochastic matrix:

$$
\mathbf{P} =
\begin{bmatrix}
\gamma^0(0) & \gamma^+(0) & 0 & 0 & \cdots & 0 \\
\gamma^{-}(\frac{1}{n}) & \gamma^0(\frac{1}{n}) & \gamma^+(\frac{1}{n}) & 0 & \cdots & 0 \\
0 & \gamma^{-}(\frac{2}{n}) & \gamma^0(\frac{2}{n}) & \gamma^+(\frac{2}{n}) & \cdots & 0 \\
0 & \cdots & \ddots & \ddots & \ddots & 0 \\
0 & \cdots & 0 & \gamma^{-}(\frac{n-1}{n}) & \gamma^0(\frac{n-1}{n}) & \gamma^+(\frac{n-1}{n}) \\
0 & \cdots & 0 & 0 & \gamma^{-}(1) & \gamma^0(1)
\end{bmatrix}
\quad (32)
$$

with $\gamma^0(c) = 1 - \gamma^+(c) - \gamma^-(c)$. This process is a random walk process. Its transition graph is strongly connected and aperiodic, hence $\mathbf{P}$ is a primitive matrix, i.e., the only absorbing class is the set of all states. This means

that in the long run, the system at time *t* can be in any
of the *n* + 1 states, and there is no stabilization [20, 21].

From Markov chain theory, the limit vector
$\pi = [\pi(0), \dots, \pi(c), \dots, \pi(1)]$ giving the probabil-
ity $\pi(c)$ to be in state $c$ in the long run is obtained as the left eigenvector of $\mathbf{P}$ associated to eigenvalue 1, i.e., $\pi$ is the solution of the system

$$ (\mathbf{P}^T - I)_z = 0 \\ 1^T z = 1 \qquad (33) $$

From now on, to avoid heavy notation, we denote $\gamma^+(k/n)$ by $\gamma^+(k)$, and similarly for $\gamma^-(k/n)$, $\pi(k/n)$, etc. We obtain

$$ \mathbf{P}^T - I = \begin{bmatrix} -\gamma^+(0) & \gamma^-(1) & 0 & 0 & 0 & \dots & 0 \\ \gamma^+(0) & -\gamma^-(1) - \gamma^+(1) & \gamma^-(2) & 0 & 0 & \dots & 0 \\ 0 & \cdots & \cdots & \cdots & 0 & \dots & 0 \\ 0 & \cdots & \gamma^+(k-1) & -\gamma^-(k) - \gamma^+(k) & \gamma^-(k+1) & \cdots & 0 \\ 0 & \cdots & 0 & \cdots & \cdots & \cdots & 0 \\ 0 & \cdots & 0 & 0 & \gamma^+(n-2) & -\gamma^-(n-1) - \gamma^+(n-1) & \gamma^-(n) \\ 0 & \cdots & 0 & 0 & 0 & \gamma^+(n-1) & -\gamma^-(n) \end{bmatrix} \quad (34) $$

Solving the system yields

$$ \begin{align*} \pi(0) &= \frac{\gamma^-(1)}{\gamma^+(0)} \pi(1) \\ \pi(1) &= \frac{\gamma^-(2)}{\gamma^+(1)} \pi(2) \\ \vdots &= \vdots \\ \pi(k) &= \frac{\gamma^-(k+1)}{\gamma^+(k)} \pi(k+1) \quad (35) \\ \vdots &= \vdots \\ \pi(n-1) &= \frac{\gamma^-(n)}{\gamma^+(n-1)} \pi(n). \end{align*} $$

This yields:

$$ \pi(k) = \frac{\gamma^+(k-1)}{\gamma^-(k)} \frac{\gamma^+(k-2)}{\gamma^-(k-1)} \cdots \frac{\gamma^+(0)}{\gamma^-(1)} \pi(0) \quad (k=1, \dots, n). \qquad (36) $$

In the case of one threshold we are able to derive the above formulas analytically. Using Eqs. (18) we obtain:

$$ \begin{align*} \pi(k) &= \frac{(1-p)(k+1)}{p(n-k)} \pi(k+1) && (k < rn-1) \\ \pi(k) &= \frac{k+1}{n-k} \pi(k+1) && (rn-1 \le k < rn) \\ \pi(k) &= \frac{p(k+1)}{(1-p)(n-k)} \pi(k+1) && (k \ge rn). \end{align*} $$

Let us find when $\pi(k)$ is increasing or decreasing. Supposing $k < rn-1$, we have:

$$ \frac{(1-p)(k+1)}{p(n-k)} \le 1 \Leftrightarrow (1-p)(k+1) \le p(n-k) \\ \Leftrightarrow k \le p(n+1)-1. $$

When $k \ge rn$, we obtain:

$$ \frac{p(k+1)}{(1-p)(n-k)} \le 1 \Leftrightarrow k \le n - p(n+1). $$

Therefore,

(i) For states below $r$, the peak is attained at

$$ \hat{c}_1 = \frac{\hat{k}_1}{n}, \text{ with } \hat{k}_1 = \lceil p(n+1) \rceil - 1. $$

Observe that when $n$ is large, this yields $\hat{c}_1 \approx p$.

(ii) For states above $r$, the peak is attained at

$$ \hat{c}_2 = \frac{\hat{k}_2}{n}, \text{ with } \hat{k}_2 = n - \lfloor p(n+1) \rfloor. $$

When $n$ is large, we obtain $\hat{c}_2 \approx 1-p$.

Depending on the relative positions of $p$ and $r$, there can be one or two peaks, as summarized as follows:

* If $r \le p$, $r \le 1-p$: peak at $\hat{c}_2$,

* if $p \le r \le 1-p$: two peaks at $\hat{c}_1, \hat{c}_2$,

* if $1-p \le r \le p$: peak at $\lceil rn \rceil / n$,

* $p \le r, 1-p \le r$: peak at $\hat{c}_1$.

In the case there are two peaks, i.e., $p \le r \le 1-p$, let us find the relative heights of the peaks. From (36), we
find, assuming $rn \notin \mathbb{N}$,

$$
\begin{align*}
\pi(\lfloor rn \rfloor) &= \pi(\hat{k}_1) \left(\frac{p}{1-p}\right)^{\lfloor rn \rfloor - \lfloor p(n+1) \rfloor + 1} \\
&\quad \times \frac{(n - \lfloor rn \rfloor + 1) \cdots (n - \lfloor p(n+1) \rfloor + 1)}{\lfloor rn \rfloor \cdots \lfloor p(n+1) \rfloor} \\
\pi(\lfloor rn \rfloor + 1) &= \pi(\hat{k}_2) \left(\frac{p}{1-p}\right)^{n - \lfloor p(n+1) \rfloor - \lfloor rn \rfloor - 1} \\
&\quad \times \frac{(\lfloor rn \rfloor + 2) \cdots (n - \lfloor p(n+1) \rfloor)}{(n - \lfloor rn \rfloor - 1) \cdots (\lfloor p(n+1) \rfloor + 1)} \\
\pi(\lfloor rn \rfloor + 1) &= \pi(\lfloor rn \rfloor) \frac{n - \lfloor rn \rfloor}{\lfloor rn \rfloor + 1}
\end{align*}
$$

Hence, assuming $p(n + 1) \notin \mathbb{N}$,

$$
\frac{\pi(\hat{k}_2)}{\pi(\hat{k}_1)} = \left(\frac{p}{1-p}\right)^{2\lfloor rn \rfloor - n + 1}. \quad (37)
$$

When $n$ is large, we obtain

$$
\frac{\pi(\hat{k}_2)}{\pi(\hat{k}_1)} \approx \left(\frac{p}{1-p}\right)^{n(2r-1)+1}. \quad (38)
$$

Observe that the peaks have equal heights when $p = 0.5$, and when $r = 0.5$, the ratio is equal to $p/(1-p)$.

Besides, we have solved numerically by SCILAB the system of equations (33), which is possible for reasonable values of $n$, and obtained its solution $\pi(k)$, $k = 0, \dots, n$. Table I shows the value of the ratio of the two peaks for various values of $p, r$ as given by Eq. (37), compared to the output of SCILAB. Fig. 5 shows the computed distribution $\pi$ for $n = 100$ for the one threshold case and also the case of the Beta distribution, compared to the histograms obtained from Monte Carlo simulations.

<table><caption>TABLE I. Example of results for $n = 10$ for different values of $p$ and $r$ under the condition $p \le r \le 1-p$. In the table are presented the theoretical ratio given by Eq. (37) (left column), as well as the values $\pi(\hat{k}_1)$ and $\pi(\hat{k}_2)$ computed numerically and the ratio between them (three rightmost columns).</caption><thead><tr><th>p</th><th>r</th><th>π(̂k<sub>2</sub>) / π(̂k<sub>1</sub>)</th><th>π(̂k<sub>1</sub>)</th><th>π(̂k<sub>2</sub>)</th><th>ratio</th></tr></thead><tbody><tr><td>0.21</td><td>0.41</td><td>3.7619048</td><td>0.064177</td><td>0.241429</td><td>3.7619052</td></tr><tr><td>0.21</td><td>0.61</td><td>0.0187835</td><td>0.2966234</td><td>0.005572</td><td>0.0187834</td></tr><tr><td>0.25</td><td>0.5</td><td>0.3333333</td><td>0.218683</td><td>0.0728942</td><td>0.3333334</td></tr></tbody></table>

We comment on these results. The Markov approach permits to obtain the stationary probability distribution of the different states, for any value of $n$, without approximation. It is found that in the long run, even if any state has a nonzero probability to be reached, some states have a much higher probability than the others to

<figure><img src="image_6.png" /><figcaption>FIG. 5. Stationary distributions of visited states for model with one threshold (upper row) with parameters (a) $r = 0.5$; $p = 0.3$, (b) $r = 0.6$; $p = 0.7$ and model with beta distribution (bottom row) with parameters (c) $\alpha = 8$; $\beta = 5$; $p = 0.1$, (d) $\alpha = 3$; $\beta = 1$; $p = 0.4$. Solid red lines are distributions obtained with Markov approach, black histograms are obtained with trajectories from Monte Carlo simulations for system of size $n = 100$ and thermalization time $t = 2 \cdot 10^6$ MCS from 100 initial conditions evenly distributed on $[0, 1]$ interval, averaged over 1000 independent runs.</figcaption></figure>

appear. In the case of one threshold, we have analytically proved the presence of one or two peaks, and their positions when $n$ is large perfectly coincides with what was predicted by the mean field approach. It is complementary to the results given by the mean-field approach, since the Markov approach is able to give the probability of occurrence of each stationary state. On the other hand, the complexity of the system of linear equations (33) induced by the Markov chain makes this approach not always tractable (e.g., with the Beta distribution). Nevertheless, we have shown that for reasonably large values of $n$ (e.g., $n = 100$), this linear system can be solved numerically, giving a perfect fit with theoretical values, as shown by Table I and with Monte Carlo simulations as well, see Fig 5.

## VII. SUMMARY AND RESEARCH DIRECTIONS FOR THE FUTURE

In this paper, we investigated the threshold model with anticonformity under asynchronous update mode, which mimics continuous time. We considered two cases: (1) homogeneous, in which all agents had the same threshold and (2) heterogeneous, in which the thresholds are given by the beta distribution function. The homogeneous case
with $r = 0.5$ is identical to the homogeneous symmetrical threshold model with anticonformity [22]. Moreover, it is almost identical to the majority-vote process [23, 24]. The only difference between the models is when the number of active and inactive agents in the neighborhood of a chosen agent is equal. In such a case, the state of the system does not change within the majority-vote model, whereas within the threshold model the change is possible. From this point of view, the threshold model with anticonformity under asynchronous updating can be treated as a generalization of a majority-vote model.

On the complete graph, the homogeneous threshold model does not give particularly interesting results. The relationship between the stationary ratio of active agents and the probability of anticonformity consists of linear dependencies, similarly as for the homogeneous symmetrical threshold model [22, 25]. The only interesting feature of this model is the discontinuity that appears at $c = r = 1 - p$. In the result, the system reaches one of two different steady states, depending on the initial conditions. Much richer behavior is observed in the heterogeneous model with thresholds given by the beta distribution function, parametrized by $\alpha, \beta$, which allows tuning the model to the homogeneous one $(\alpha, \beta \to \infty)$ to maximally heterogeneous (i.e. described by the uniform distribution function).

A particularly interesting behavior is obtained if at least one of the shape parameters $\alpha$ or $\beta$ is larger than one and both parameters are finite. In this case PDF has a shape that reminds those of actual trait manifestation in behavior, as reported by psychologists [17], i.e., uni modal, not necessarily symmetrical, function with maximum at the value $0 < r < 1$. In such a case a phase transition appears, which is continuous for $\alpha = \beta$, and discontinuous otherwise. In the latter case, the transition involves phenomena typical of social systems, such as social hysteresis [26] and the critical mass [27].

The future research on the model can be conducted in several directions, related to the following questions:

* How the results would change if the threshold for anticonformity would be different than that for conformity? This question is inspired by the work on the *q*-voter model with generalized anticonformity

[28]. In the *q*-voter model such a generalization resulted in switching from continuous to discontinuous phase transitions for some values of parameters. The question is if the same phenomena would be observed for the threshold model.

* How the structure of a network would influence the results? This question is inspired by the work on the symmetrical threshold [25]. It was shown that on random graphs with the degree observed empirically for social networks, the largest social hysteresis is observed for $r \in (0.65, 0.85)$. This was a meaningful result from the social point of view and thus it would be desirable to check if it appears also in the asymmetric model studied here.

* How the results would change if the quenched approach to anticonformity would be used. In this version of the model, we used the annealed approach, in the sense that each agent could anticonform (with probability $p$) or conform (with probability $1-p$). However, we could use also the quenched approach, in which a fraction $p$ of agents are permanently anticonformists. This question is inspired by the work on the *q*-voter model with nonconformity under quenched and annealed approaches [29]. It was shown that on the complete graph both approaches give the same result for the *q*-voter model with anticonformity, whereas different for the model with independence. The question is to what extend this result is universal.

## ACKNOWLEDGMENTS

This work was supported by the National Science Center (NCN, Poland) through Grant No. 2019/35/B/HS6/02530

## ACKNOWLEDGMENTS

This research was supported by the National Science Center (NCN, Poland) grant number 2019/35/B/HS6/02530.

[1] J. P. Gleeson, Binary-state dynamics on complex networks: Pair approximation and beyond, Phys. Rev. X 3, 021004 (2013).

[2] M. Granovetter, Threshold models of collective behavior, American Journal of Sociology 83, 1420-1443 (1978).

[3] D. Watts, A simple model of global cascades on random networks, Proceedings of the National Academy of Sciences of the United States of America 99, 5766 (2002).

[4] V. V. Breer, Models of tolerant threshold behavior (from T. Schelling to M. Granovetter), Automation and Remote Control 78, 10.1134/S0005117917070128 (2017).

[5] A. Jędrzejewski and K. Sznajd-Weron, Statistical Physics Of Opinion Formation: Is it a SPOOF?, C. R. Physique 20, 244 (2019).

[6] T. C. Schelling, Micromotives and macrobehavior (New York: Norton, 1978).

[7] P. Dodds and D. Watts, Universal behavior in a generalized model of contagion, Physical Review Letters 92, 218701 (2004).

[8] J. S. Juul and M. A. Porter, Hipsters on networks: How a minority group of individuals can lead to an antiestablishment majority, Phys. Rev. E 99, 022313 (2019).
[9] E. Lee and P. Holme, Social contagion with degree-dependent thresholds, Physical Review E 96, 10.1103/PhysRevE.96.012315 (2017).

[10] M. Grabisch and F. Li, Anti-conformism in the threshold model of collective behavior, Dynamic Games and Applications 10, 10.1007/s13235-019-00332-0 (2020).

[11] C. Castellano, S. Fortunato, and V. Loreto, Statistical physics of social dynamics, Rev. Mod. Phys. 81, 591 (2009).

[12] S. Galam and A. Martins, Two-dimensional Ising transition through a technique from two-state opinion-dynamics models, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 91, 10.1103/PhysRevE.91.012108 (2015).

[13] T. Raducha, M. Wilinski, T. Gubiec, and H. Stanley, Statistical mechanics of a coevolving spin system, Physical Review E 98, 10.1103/PhysRevE.98.030301 (2018).

[14] M. Calvelli, N. Crokidakis, and T. J. Penna, Phase transitions and universality in the sznajd model with anticonformity, Physica A 513, 518 (2019).

[15] A. R. Vieira, A. F. Peralta, R. Toral, M. S. Miguel, and C. Anteneodo, Pair approximation for the noisy threshold q-voter model, Phys. Rev. E 101, 052131 (2020).

[16] S. H. Strogatz, *Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (2nd ed.)* (CRC Press, 2015).

[17] W. Fleeson and P. Gallagher, The implications of big five standing for the distribution of trait manifestation in behavior: Fifteen experience-sampling studies and a meta-analysis, Journal of Personality and Social Psychology 97, 1097 (2009).

[18] V. Yang, M. Galesic, H. McGuinness, and A. Harutyunyan, Dynamical system model predicts when social learn-

ers impair collective performance, Proceedings of the National Academy of Sciences of the United States of America 118, 10.1073/pnas.2106292118 (2021).

[19] S. Banisch, R. Lima, and T. Araújo, Agent based models and opinion dynamics as markov chains, Social Networks 34, 549 (2012).

[20] J. G. Kemeny and J. L. Snell, *Finite Markov chains* (Springer Verlag, 1976).

[21] E. Seneta, *Non-negative Matrices and Markov Chains* (Springer, 2006).

[22] B. Nowak and K. Sznajd-Weron, Homogeneous symmetrical threshold model with nonconformity: Independence versus anticonformity, Complexity 2019, 1 (2019).

[23] T. M. Liggett, *Interacting Particle Systems* (Springer, 1985).

[24] M. de Oliveira, Isotropic majority-vote model on a square lattice, Journal of Statistical Physics 66, 273 (1992).

[25] B. Nowak and K. Sznajd-Weron, Promoting discontinuous phase transitions by the quenched disorder within the multistate q-voter model (2021), arXiv:2106.11238 [physics.soc-ph].

[26] M. Scheffer, F. Westley, and W. Brock, Slow response of societies to new problems: Causes and costs, Ecosystems 6, 493 (2003).

[27] D. Centola, J. Becker, D. Brackbill, and A. Baronchelli, Experimental evidence for tipping points in social convention, Science 360, 1116 (2018).

[28] A. Abramiuk-Szurlej, A. Lipiecki, J. Pawłowski, and K. Sznajd-Weron, Discontinuous phase transitions in the q-voter model with generalized anticonformity on random graphs, Scientific Reports 11, 1 (2021).

[29] A. Jędrzejewski and K. Sznajd-Weron, Person-situation debate revisited: Phase transitions with quenched and annealed disorders, Entropy 19, 415 (2017).