Collective dynamics of random Janus oscillator networks

TL;DR

This paper investigates the collective behavior of Janus oscillators on various large networks, revealing persistent complex dynamics, coexistence of states, and abrupt transitions, with unique features in scale-free networks.

Contribution

It demonstrates that Janus oscillator dynamics persist in large random networks and uncovers novel phenomena like state coexistence and abrupt transitions across different topologies.

Findings

Rich dynamics persist in large networks.

Coexistence of synchronized and oscillatory states.

Abrupt transitions in order parameters.

Abstract

Janus oscillators have been recently introduced as a remarkably simple phase oscillator model that exhibits non-trivial dynamical patterns -- such as chimeras, explosive transitions, and asymmetry-induced synchronization -- that once were only observed in specifically tailored models. Here we study ensembles of Janus oscillators coupled on large homogeneous and heterogeneous networks. By virtue of the Ott-Antonsen reduction scheme, we find that the rich dynamics of Janus oscillators persists in the thermodynamic limit of random regular, Erd\H{o}s-R\'enyi and scale-free random networks. We uncover for all these networks the coexistence between partially synchronized state and a multitude of states displaying global oscillations. Furthermore, abrupt transitions of the global and local order parameters are observed for all topologies considered. Interestingly, only for scale-free networks, it is found that states displaying global oscillations vanish in the thermodynamic limit.

Equations

$${\dot{\theta }}_i={\omega }_i+\sum _{j=1}^{2N}{W}_{ij}\sin ({\theta }_j-{\theta }_i),(i=1,...,2N)$$

$$\mathbf{W}=\left[\begin{array}[]{cc}0&\beta\mathbf{\mathbf{I}+\sigma A}\\ \beta\mathbf{I}+\sigma\mathbf{A}^{{\rm T}}&0\end{array}\right].$$

$$R_i=\sum _{j=1}^{2N}{W}_{ij}{e}^{\mathrm{i}{\theta }_j},$$

$${\dot{\theta }}_i={\omega }_i+\text{Im}[e^{-\mathrm{i}{\theta }_i}{R}_i].$$

$$\frac{\partial {\rho }_i}{\partial t}+\frac{\partial }{\partial {\theta }_i}({\rho }_i{\dot{\theta }}_i)=0.$$

$${\rho }_i({\theta }_i,{\omega }_i,t)=\frac{g({\omega }_i)}{2\pi }\left[1+\sum _{n=1}^{\infty }{\hat{\alpha }}_i^n({\omega }_i,t)e^{\text{i}n{\theta }_i}+c.c\right]$$

$$\frac{d{\hat{\alpha }}_i}{dt}+\mathrm{i}{\hat{\alpha }}_i{\omega }_i+\frac{1}{2}\left[{\hat{\alpha }}_i^2{R}_i-R_i^{\ast }\right]=0(i=1,...2N),$$

$$R_i=\sum _{j=1}^{2N}{W}_{ij}{\int }_{-\infty }^{\infty }{\int }_0^{2\pi }{\rho }_j({\theta }_j,{\omega }_j,t)e^{\text{i}{\theta }_j}d{\theta }_{j}d{\omega }_j.$$

$$R_i=\sum _{j=1}^{2N}{W}_{ij}{\int }_{-\infty }^{\infty }{\hat{\alpha }}_j^{\ast }({\omega }_j,t)g({\omega }_j)d{\omega }_j.$$

$${\hat{\alpha }}_1={\hat{\alpha }}_2=\cdots ={\hat{\alpha }}_N\equiv {\alpha }_1$$

$${\hat{\alpha }}_{N+1}=\cdots ={\hat{\alpha }}_{2N}\equiv {\alpha }_2.$$

$$R_i=\{\begin{array}{ll}{\alpha }_2^{\ast }(\beta +k\sigma )\equiv R_1& \text{ if }i=1,...,N;\\ {\alpha }_1^{\ast }(\beta +k\sigma )\equiv R_2& \text{ if }i=N+1,...,2N.\end{array}$$

$${\dot{\alpha }}_{1,2}+\mathrm{i}\omega {\alpha }_{1,2}+\frac{1}{2}\left[{\alpha }_{1,2}^2{R}_{1,2}-R_{1,2}^{\ast }\right]=0,$$

$${\dot{r}}_1$$

$${\dot{r}}_2$$

$$\dot{\delta }$$

$${\sigma }_{c_1}=-\frac{\Delta /2+\beta }{k}\text{ and }{\sigma }_{c_2}=\frac{\Delta /2-\beta }{k}.$$

$$R=\frac{1}{\sqrt{2}}\sqrt{1\pm \sqrt{1-\frac{{\Delta }^2}{4(\beta +k\sigma {)}^2}}},$$

$${\sigma }_{c_3}\equiv -\frac{(\Delta +\beta )}{k}<\sigma <\frac{\Delta -\beta }{k}\equiv {\sigma }_{c_4}.$$

$${\dot{b}}_{\perp }=\gamma b_{\perp }(b_{\perp }-b_{\parallel }-1)\cos \delta .$$

$${\lambda }_{\perp }=-(\beta +k\sigma )(1+b_0)\cos {\delta }_0.$$

$${\dot{\alpha }}_{k,1}+\mathrm{i}\omega {\alpha }_{k,1}+\frac{1}{2}[{\alpha }_{k,1}^2\left(\beta {\alpha }_{k,2}^{\ast }+\sigma \frac{k}{⟨k⟩}\sum _{k^{\prime }}{k}^{\prime }{p}_{k^{\prime }}{\alpha }_{k^{\prime },2}^{\ast }\right)$$

$$-\left(\beta {\alpha }_{k,2}+\sigma \frac{k}{⟨k⟩}\sum _{k^{\prime }}{k}^{\prime }{p}_{k^{\prime }}{\alpha }_{k^{\prime },2}\right)]$$

$$\dot{\delta \alpha }=-\mathrm{i}\omega \delta \bar{\alpha }+\frac{1}{2}\left[\beta +\sigma \frac{⟨k^2⟩}{⟨k⟩}\right]\delta \alpha$$

$$\dot{\delta \bar{\alpha }}=-\mathrm{i}\omega \delta \alpha -\frac{1}{2}\left[\beta +\sigma \frac{⟨k^2⟩}{⟨k⟩}\right]\delta \bar{\alpha },$$

$${\sigma }_{c_1}^{(\mathrm{net})}\equiv -(\Delta +\beta )\frac{⟨k⟩}{⟨k^2⟩}<\sigma <\frac{⟨k⟩}{⟨k^2⟩}(\Delta -\beta )\equiv {\sigma }_{c_2}^{(\mathrm{net})}.$$

Full text

Research on coupled oscillators in the past decade has been marked by the discovery of many intriguing patterns in the collective behavior of networks Arenas et al. (2008); Rodrigues et al. (2016). Notable examples of such patterns are chimeras Panaggio and Abrams (2015), states in which populations of synchronous are asynchronous oscillators coexist; explosive synchronization transitions Gómez-Gardenes et al. (2011); Rodrigues et al. (2016), which appear as a consequent of constraints in the natural frequency assignment; and asymmetry-induced synchronization Nishikawa and Motter (2016); *zhang2017asymmetry, a state in which synchrony is counter-intuitively favored by oscillator heterogeneity. In all these cases, phase oscillator models had to be specially designed so that those non-trivial states could be scrutinized. Very recently, however, Nicolaou et al. Nicolaou et al. (2019) defined an oscillator model coined as Janus oscillators; the name is inspired in the homonym two-faced god of Roman mythology and reflects the two-dimensional character of an isolated oscillator – each “face” of a Janus unit consists of a Kuramoto oscillator, whose natural frequency has the same absolute value but opposite sign to the frequency of its counter-face. When coupled on one-dimensional regular graphs, Janus oscillators have been found to exhibit a striking rich dynamical behavior that encompasses the co-occurrence of several dynamical patterns, in spite of the simplicity of the topology and the oscillator model itself Nicolaou et al. (2019).

The Janus model was introduced as a potential model for biological systems such as the Chlamydomonas cells with counterrotating flagella Friedrich and Jülicher (2012); Wan and Goldstein (2016). It is thus important to understand the dynamics of a Janus system on related (realistic) topologies. Here we pose the question of whether the observed 1D rich collective dynamics exists on more complex networks of Janus oscillators. To address this issue, we employ the Ott-Antonsen (OA) ansatz Ott and Antonsen (2008) and obtain a reduced set of equations describing the system’s evolution. From this reduced representation we find that, indeed, peculiar patterns of synchrony persist when Janus oscillators are placed on random regular, Erdős-Rényi (ER) and scale-free (SF) random networks. We provide analytical and numerical evidence that the multitude of states in Janus dynamics is a consequence of the coexistence of infinite neutrally stable limit-cycle trajectories, which we denominate “breathing standing-waves”. Co-occurrence between classical partially synchronized and standing-waves are also reported. We further show that for high average degrees the collective states of ER networks are accurately described by the reduced system obtained for random regular ones. Interestingly, we demonstrate that the coupling range in which global oscillations are possible vanishes in the thermodynamic limit of SF networks.

We begin by defining the dynamics of $N$ Janus oscillators Nicolaou et al. (2019) on heterogeneous networks as

[TABLE]

where the $2N\times 2N$ matrix $\mathbf{W}$ is defined as

[TABLE]

Natural frequencies are assigned as $\omega_{i}=\omega_{0}+\Delta/2$, for $i=1,...,N$; and $\omega_{i}=\omega_{0}-\Delta/2$, for $i=N+1,...,2N$, where $\Delta$ is the frequency mismatch and $\omega_{0}$ is the average frequency, which we assume $\omega_{0}=0$. System (1) is analogous to a bipartite network or a multilayer network in which oscillators belonging to the same group do not interact with one another, while connections between groups are encoded in matrix $\mathbf{A}$. Notice also that interactions between oscillators are weighted by the coupling strength $\sigma$, except for oscillators with indexes $(i,i+N)$, $i\in[1,N]$ – these pairs of nodes interact with coupling strength $\beta$.

By defining the local order parameters

[TABLE]

Eqs. 1 are then decoupled as

[TABLE]

Following Barlev et al. (2011), we consider an ensemble of systems defined by Eq. 1 with fixed coupling matrix $\mathbf{W}$. In this formulation, we describe the system (1) at a given time step $t$ by the joint probability density $\rho_{2N}(\bm{\theta},\bm{\omega},t)$, where $\bm{\theta}=(\theta_{1},...,\theta_{2N})$ is the vector containing the phases at time $t$, and $\bm{\omega}=(\omega_{1},...,\omega_{2N})$ is the time-independent vector with the natural frequencies of the individual oscillators. The evolution of the join probability $\rho_{2N}$ is then dictated by $\partial_{t}\rho_{2N}+\sum_{i=1}^{2N}\partial_{\theta_{i}}(\rho_{2N}\dot{\theta}_{i})=0$, where $\dot{\theta}_{i}$ is given by Eq. 4. Let us suppose that frequencies $\omega_{j}$ are distributed according to a generic function $g(\omega_{j})$. Multiplying the continuity equation of $\rho_{2N}$ by $\Pi_{j\neq i}d\omega_{j}d\theta_{j}$ and integrating, we obtain the evolution equation for the marginal oscillator density $\rho_{i}(\theta_{i},\omega_{i},t)=\int\int\rho_{2N}\Pi_{j\neq i}d\omega_{j}\theta_{j}$; that is

[TABLE]

By expanding $\rho_{i}$ in Fourier series and applying the OA ansatz to its coefficients we have Ott and Antonsen (2008); Barlev et al. (2011)

[TABLE]

Inserting the previous equation into Eq. 5, we obtain the evolution for coefficients $\hat{\alpha}_{i}$:

[TABLE]

where, in this ensemble approach, the coefficients $R_{i}$ are calculated as

[TABLE]

Inserting Eq. 6 in the previous equation yields

[TABLE]

Let us now consider the case in which each subpopulation of the Janus coupling scheme is a random regular network with degree $k$. More precisely, each oscillator of the subpopulation that rotates with frequency $\omega_{1}=\Delta/2$ is randomly connected to $k$ oscillators of subpopulation 2 ($\omega_{2}=-\Delta/2$), and vice-versa. In this scenario, since oscillators within each group are identical – i.e. $g(\omega_{i})=\delta(\omega_{i}-\Delta/2)$, for $i=1,...,N$; and $g(\omega_{i})=\delta(\omega_{i}+\Delta/2)$, for $i=N+1,...,2N$–, we assume the following solution for coefficients $\hat{\alpha}_{i}$:

[TABLE]

Hence, the local order parameters $R_{i}$ are reduced to:

[TABLE]

Inserting solutions (10) and (11) in Eq. 7, we obtain the reduced set of equations

[TABLE]

which in terms of the coordinates $\alpha_{1,2}=r_{1,2}e^{i\psi_{1,2}}$ are written as

[TABLE]

where $\delta=\psi_{1}-\psi_{2}$ is the phase-lag between subpopulations. Variables $r_{1}$ and $r_{2}$ turn out to be the order parameters measuring the level of synchronization within each subpopulation in the Janus system; the traditional Kuramoto order parameter evaluating the global synchrony is obtained through $R(t)=\frac{1}{2}|r_{1}e^{i\psi_{t}(t)}+r_{2}e^{i\psi_{2}(t)}|$. States that emerge from system (13) are summarized as: (1) a partially synchronized state in which $r_{1,2}=1$, while the subpopulations remain separated by a constant phase-lag $\delta$ (hence, $R<1$); (2) a standing-wave (SW) state, where the bulks of the two fully synchronized populations ($r_{1,2}=1$) rotate in opposite directions yielding a incessantly rotating $\delta$; (3) a distinct form of SW emerges when $0<r_{1,2}<1$: along with the increment or decrease in phase-lag $\delta$, the order parameters $r_{1,2}$ exhibit a breathing behavior, as depicted in the simulation shown Fig. 1. Henceforth we refer to this state as “breathing SW”. As we shall see, the classical incoherent state remains unstable for all coupling values.

In order to uncover the conditions for the existence of the partially synchronized state, we set $r_{1}(t)=r_{2}(t)\equiv r(t)$ in Eqs. 13, leading to $\dot{r}=\frac{\gamma r}{2}(1-r^{2})\cos\delta$ and $\dot{\delta}=-\Delta-\gamma(1+r^{2})\sin\delta$. Imposing $\dot{r}=0$, we notice that $r=1$ is always a fixed point. Inserting the latter solution in the equation for $\dot{\delta}=0$, we find that $\sin\delta=-\Delta/2\gamma$. Thus, the partially synchronized regime remains stable when $-\Delta/2|\gamma|\leq 1$ is satisfied. In terms of coupling $\sigma$, we then write these critical conditions as

[TABLE]

Couplings $\sigma_{c_{1}}$ and $\sigma_{c_{2}}$ determine the coupling range where the partially synchronized state exists. The total order parameter $R$ is then given by

[TABLE]

where the “-” branch is stable for $\sigma\leq\sigma_{c_{1}}$, whereas the “+” branch is stable in the region $\sigma\geq\sigma_{c_{2}}$. For $\sigma_{c_{1}}<\sigma<\sigma_{c_{2}}$, the limit cycle solution of $\dot{\delta}$ holds and SW states emerge with perfectly synchronized subpopulations ($r_{1,2}=1$).

A linear stability analysis of the incoherent state ($\alpha_{1,2}=0$) in Eqs. 12 reveals that the eigenvalues of the Jacobian matrix become purely imaginary at $|\beta+k\sigma|=\Delta$. Therefore, limit cycle solutions arise for

[TABLE]

Figure 2(a) outlines the critical conditions given by Eqs. 14 and  16 in the plane spanned by couplings $\beta$ and $\sigma$. As it can be seen, the partially synchronized state is favored by extreme values of both couplings, whereas states with oscillating synchrony appear for intermediate values in the parameter space. In Figs. 2 (b) and (c) we visualize the evolution of the local and global order parameters over two vertical sections of the diagram in (a), namely for $\beta=0.25$ and $-0.3$, respectively. Supposing we initiate the system with a negative $\sigma$ in the “Partial sync” region ($\sigma<\sigma_{c_{4}}$), the total order parameter $R$ collapses in the solution given by Eq. 15. As $\sigma$ is further increased, at $\sigma=\sigma_{c_{2}}$ the unstable and stable branches of $R$ merge via a saddle-node infinite-period bifurcation, whereby the limit cycle solution of the SW arises (see Figs. 2 (b) and (c)). Upon further continuation of $\sigma$, a saddle-node appears again at $\sigma=\sigma_{c_{1}}$ and the system is brought back to the partially synchronized state.

Besides the branch of $R$ obtained under the symmetry condition $r_{1}=r_{2}=1$, the numerical evolution of the reduced system predicts the existence of several other curves, which are upper bounded by $R$ of the $r_{1,2}=1$ solution for $\sigma\in[\sigma_{c_{4}},\sigma_{c_{3}}]$. Insights about the nature of such states can be gained by investigating the stability of the SW state under perturbations transversal to the symmetric manifold $r_{1}=r_{2}$ (Martens et al., 2009). By defining the transversal and longitudinal coordinates $r_{\perp}=(r_{1}-r_{2})/2$ and $r_{\parallel}=(r_{1}+r_{2})/2$, we have $\dot{r}_{\perp}=\frac{\gamma}{2}(r_{\perp}^{2}-r_{\parallel}^{2}-1)r_{\perp}\cos\delta$, which in terms of variables $b_{\parallel}=r_{\parallel}^{2}$ and $b_{\perp}=r_{\perp}^{2}$ reads

[TABLE]

Linearization at a point $(b_{0},\delta_{0})$ lying on a limit cycle solution of the manifold $r_{1}=r_{2}$ yields the variation equation $\dot{\delta b}_{\perp}=\lambda_{\perp}\delta b_{\perp}$, where

[TABLE]

By averaging the previous equation over a period of oscillation and using the periodicity of the limit cycle ($\langle(d/dt)\ln b_{0}\rangle=0$) we find $\langle\lambda_{\perp}\rangle=-2(\beta+k\sigma)\langle\cos\delta_{0}\rangle$. Numerical calculations with original system (1) for extensive parameter combinations show that $\langle\cos\delta_{0}\rangle\approx 0$ (and, consequently, $\langle\lambda_{\perp}\rangle\approx 0$) in the region $\sigma\in[\sigma_{c_{2}},\sigma_{c_{1}}]$, suggesting then that the limit cycle solution of the SW is neutrally stable. Although our numerical estimate does not give us an exact proof of the stability of the limit cycle solution, it sheds light on the existence of the multitude of curves observed in Fig. 2(b) and (c). Essentially, any perturbation of the SW state leads to a new limit cycle with $r_{1,2}<1$, since nearby trajectories are not attracted nor repelled, explaining the origin of the numerous solutions encountered in the coupling region encompassed by $\sigma_{c_{4}}$ and $\sigma_{c_{3}}$ in Fig. 2. Notice also that the lower branches in the region $\sigma\in[\sigma_{c_{4}},\sigma_{c_{3}}]$ do not correspond precisely to the classical incoherent state, but rather represent limit cycles solutions with small amplitudes $r_{1,2}$. It is noteworthy mentioning also that although we have considered negative and positive couplings, we have not observed in the populations of Janus oscillators states akin to traveling waves and $\pi$-states, which are collective phenomena that are characteristic of the interplay between attractive and repulsive interactions Hong and Strogatz (2011); Sonnenschein et al. (2015).

The theory developed for random regular networks can also provide insights on the dynamics of networks with mildly heterogeneous degree distributions. In Fig. 3(a) we superimpose numerical results for ER networks with the branches for the partially synchronized states (Eq. 15) and critical conditions given by Eq. 14 and 16. Interestingly, we see that the dependence of the order parameters is reproduced with good precision with the expressions derived for simpler networks. Boundaries enclosing the breathing SW states in the random regular network also delineate the region with global oscillations for the ER network. Notice also that $\sigma_{c_{1,2}}$ again mark the limits of the partially synchronize branch of $R$; however, no state analogous to the perfectly symmetric SW ($r_{1,2}=1$) is observed in $\sigma_{c_{2}}<\sigma<\sigma_{c_{1}}$ for ER networks.

Let us take a step further in the analysis of heterogeneous structures and consider general uncorrelated networks with degree distribution $p_{k}$. In this case, we assume that nodes with same degree $k$ admit the same solution, i.e., $\alpha_{i}=\alpha_{k}$, if $k_{i}=k$. Thus, Eqs. 7 are reduced to

[TABLE]

where $\alpha_{k,1}$ describes the dynamics of oscillators with degree $k$ and frequency $\omega_{i}=\Delta/2$. Equations for coefficients $\alpha_{k,2}$ standing for the second face of Janus oscillators ($\omega_{i}=-\Delta/2$) are obtained accordingly. Linearizing the system around $\alpha_{k,1}=\delta\alpha_{k,1}\ll 1$ yields the following variational system

[TABLE]

where $\delta\alpha$ is a small perturbation of the complex order parameter $\alpha=\frac{1}{2\langle k\rangle}\sum_{k}kp_{k}(\alpha_{k,1}+\alpha_{k,2})$, and $\delta\bar{\alpha}$ is the analogous quantity for the parameter measuring the difference in the internal synchrony in Janus oscillators, i.e., $\bar{\alpha}=\frac{1}{2\langle k\rangle}\sum_{k}kp_{k}(\alpha_{k,1}-\alpha_{k,2})$. The eigenvalues of the Jacobian matrix of system (19) become purely imaginary at $\Delta=|\beta+\sigma\langle k^{2}\rangle/\langle k\rangle|$. Therefore, we predict the appearance of states with global oscillations in the range

[TABLE]

We check the predictions of the equation above in Fig. 3(b) for SF networks with power-law exponent $\gamma=2.25$. At first sight, it seems that the condition $\sigma_{c_{1}}^{(\rm{net})}$ provides an inaccurate estimation of the region where the order parameters $R$ and $r_{1,2}$ are expected to exhibit an erratic behavior, suggesting perhaps that finite-size effects could be behind the deviation between $\sigma_{c_{1}}^{(\rm net)}$ and the point $\sigma\simeq 0.17$ at which the branch of $R$ collapses to values $R\approx 0$. However, Eq. 20 refers to the coupling range in which multiple oscillating states are expected to emerge. Visualizing in Fig. 3(c) the evolution of the mean-field frequency $\dot{\Theta}$, we observe that $\sigma_{c_{1,2}}^{(\rm{net})}$ actually define very accurately the boundaries of the states with multiple oscillating solutions ($\sigma_{c_{1}}^{\rm{(net)}}<\sigma<\sigma_{c_{2}}^{\rm{(net)}}$). Given the dependence on $\langle k\rangle/\langle k^{2}\rangle$, one further envisions from Eq. 20 the absence of such oscillating states in the thermodynamic limit for SF networks with $2<\gamma\leq 3$, since $\sigma_{c_{1,2}}^{\rm{(net)}}$ are expected to vanish as $N\rightarrow\infty$, similarly to the classical coupling for the onset of synchronization in such structures Peron et al. (2019); Rodrigues et al. (2016).

In conclusion, we have explored the collective dynamics of Janus oscillators on large homogeneous and heterogeneous random networks. By employing the scheme provided by the OA ansatz, we obtained, for random regular networks, a reduced set of equations whereby critical points of the dynamics were revealed. We found that several collective behaviors coexist for intermediate coupling values, elucidating the findings in Nicolaou et al. (2019). Although initially obtained for homogeneous networks, we verified that the solutions of the reduced system fitted accurately numerical experiments for dense ER networks. By analyzing the stability of general uncorrelated networks, we further verified that the coupling range for which global oscillations are possible shrinks in the thermodynamic limit of SF networks. It is pertinent to mention that the accuracy of the OA ansatz in predicting the transition points is deteriorated for $\sigma$ and $\beta$ values beyond the region depicted in Fig. 2. Deviations from the temporal signature yielded by the reduced system for solutions of $r_{1,2}$ were also observed in simulations for some couplings $(\beta,\sigma)$ in the breathing SW area.

All in all, we provided the first theoretical and numerical analysis of ensembles of Janus oscillators on homogeneous and heterogeneous random networks. As such, our work raises further interesting questions about the study initiated by Nicolaou et al. Nicolaou et al. (2019). For instance, future investigations should target the dynamics on sparse and correlated networks – situations in which the ensemble approach in Barlev et al. (2011) and mean-field techniques are inaccurate in predicting tipping points of the system – as well as limitations of the OA manifold in capturing the Janus dynamics.

TP acknowledges FAPESP (Grants No. 2016/23827-6 and 2018/15589-3). DE acknowledges Kadir Has University internal Scientific Research Grant (BAF). FAR acknowledges CNPq (Grant No. 305940/2010-4) and FAPESP (Grants No. 2016/25682-5 and grants 2013/07375-0) for the financial support given to this research. YM acknowledges partial support from the Government of Aragon, Spain through grant E36-17R (FENOL), by MINECO and FEDER funds (FIS2017-87519-P) and by Intesa Sanpaolo Innovation Center. This research was carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (grant 2013/07375-0). The funders had no role in study design, data collection and analysis, or preparation of the manuscript.