Non-ideal rheology of semidilute bacterial suspensions

TL;DR

This paper investigates the complex rheological behavior of semidilute bacterial suspensions using kinetic theory, revealing non-linear effects, phase transitions, and shear stress oscillations that challenge idealized models.

Contribution

It introduces models considering binary interactions with steric and short-range forces, demonstrating non-ideal rheology and phase transitions in bacterial suspensions.

Findings

Ordered phase develops with increasing density.

Shear stress oscillates with large amplitudes at certain shear rates.

Non-linear effects can increase shear stress in pusher swimmers.

Abstract

The rheology of semidilute bacterial suspensions is studied with the tools of kinetic theory, considering binary interactions, going beyond the ideal gas approximation. Two models for the interactions are considered, which encompass both the steric and short range interactions. In these, swimmers can either align polarly regardless of the state previous to the collision or they can align axially, being possible the end up antiparallel if the relative angle between directors is large. In both cases, it is found that an ordered phase develops when increasing the density, where the shear stress oscillates with large amplitudes, when a constant shear rate is imposed. This oscillation disappears for large shear rates in a continuous or discontinuous transition, depending if the aligning is polar or axial, respectively. For pusher swimmers these non-linear effects can produce an increase on the shear stress, contrary to the prediction of viscosity reduction made for the dilute regime with the ideal gas approximation.

Figures

Equations

$$\dot{x}$$

$$\dot{\theta }$$

$$\frac{\partial \Psi }{\partial t}+\frac{\partial (\dot{x}\Psi )}{\partial x}+\frac{\partial (\dot{y}\Psi )}{\partial y}+\frac{\partial (\dot{\theta }\Psi )}{\partial \theta }=D_{\text{r}}\frac{{\partial }^2\Psi }{\partial {\theta }^2}+J,$$

$$J=g{\int }_{-\pi }^{\pi }dw\left[\Psi (\theta +w/2)\Psi (\theta -w/2)-\Psi (\theta )\Psi (\theta -w)\right],$$

$${\Sigma }_{xy}^{\text{a}}={\sigma }_0\int \Psi \cos \theta \sin \theta \mathrm{d}\theta ,$$

$$\frac{d{b}_n}{ds}=$$

$$+2\pi \sum _l{b}_l{b}_{n-l}\text{sinc}\left(n\pi /2-l\pi \right),$$

$$⟨\hat{\mathbf{p}}⟩$$

$$\frac{d{b}_1}{ds}$$

$$\frac{d{b}_2}{ds}$$

$$b_1(s_0,s_1)=A{e}^{i{\omega }_0{s}_0+i{\omega }_1{s}_1}+B{e}^{-i{\omega }_0{s}_0-i{\omega }_1{s}_1},$$

$$b_2(s_0,s_1)=F{e}^{2i{\omega }_0{s}_0+2i{\omega }_1{s}_1}+G{e}^{-2i{\omega }_0{s}_0-2i{\omega }_1{s}_1}+H,$$

$$H_0=-\frac{i\beta \zeta {\rho }_c}{4\pi }\frac{1}{{\rho }_c+4-i\zeta }.$$

$$\left(\begin{array}{cc}-6i(\zeta -2{\omega }_0)& 32H_0+3i\beta \zeta \\ 32H_0^{\ast }-3i\beta \zeta & 6i(\zeta +2{\omega }_0)\end{array}\right)\left(\begin{array}{c}A\\ B^{\ast }\end{array}\right)=0,$$

$$B^{\ast }=\frac{6i(\zeta -2{\omega }_0)}{32H_0+3i\beta \zeta }A\equiv lA$$

$${\omega }_0^2=\frac{|32H_0+3i\beta \zeta {|}^2-36{\zeta }^2}{4}.$$

$$F$$

$$G$$

$$H_1$$

$$\equiv h_3\frac{AB}{ϵ}+h_4.$$

$$ϵ\frac{dA}{d{s}_1}=A\left[ϵ\left(1-i{\omega }_1-\frac{8}{3}{h}_{4}l\right)-\frac{8}{3}\left(h_1+h_3|l{|}^2\right)|A{|}^2\right],$$

$${\Sigma }_{xy}^{\text{a}}$$

$${\Sigma }_{xy}^{\text{a,osc}}$$

$${\Sigma }_{xy}^{\text{a,mean}}$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!+\left[48+32\epsilon+\zeta^{2}(1+2\epsilon)\right]\rho_{c}+4(3+\epsilon)\rho_{c}^{2}+\rho_{c}^{3}\bigg{\}}.$$

$$J=g\bigg{\{}\int_{-\pi/2}^{\pi/2}dw\big{[}\Psi(\theta+w/2)\Psi(\theta-w/2)\\ +\Psi(\theta+w/2)\Psi(\theta+\pi-w/2)\big{]}-\int_{-\pi}^{\pi}\Psi(\theta)\Psi(\theta-w)\bigg{\}}.$$

$$\frac{d{b}_n}{ds}=$$

$$+\pi \sum _l{b}_l{b}_{n-l}\text{sinc}\left(n\pi /4-l\pi /2\right)\left[1+(-1{)}^{n-l}\right],$$

Full text

I Introduction

Bacterial suspensions have attracted great attention in recent years as one of the quintessential examples of active matter Berg (2008); Vicsek and Zafeiris (2012); Marchetti et al. (2013). They exhibit striking features which are absent in their passive counterparts, such as collective oscillatory behavior Chen et al. (2017), active turbulence Wensink et al. (2012); Cisneros et al. (2007), hydrodynamic instabilities Saintillan and Shelley (2008); Pahlavan and Saintillan (2011), and exotic rheological properties Gachelin et al. (2013); López et al. (2015); Sokolov and Aranson (2009); Rafaï et al. (2010); Saintillan (2010) (for a recent review on the rheology, see Ref. Saintillan (2018)). Mechanically, swimmers can be classified according to the sign of the force dipole they exert on the fluid as pushers for extensile dipoles and pullers for contractile ones. Experiments on Escherichia coli suspensions (pushers) revealed a decrease in the total viscosity with the swimmer concentration Gachelin et al. (2013); Sokolov and Aranson (2009), reaching even “superfluid” states at moderate concentrations López et al. (2015), whereas experiments on Chlamydomonas reinhardtii (pullers) exhibited a net increase on the viscosity Rafaï et al. (2010). These modifications on the suspension viscosity have been accounted for through kinetic theories Saintillan (2010) in the dilute regime, i.e. when all the swimmer–swimmer interactions are neglected in the ideal gas rheology. It is found that for small shear rates, the active contribution to the viscosity is proportional to the dipole strength and bacterial concentration, being negative for pushers and positive for pullers in accordance with the experiments. This effect depends on the swimmer shape, vanishing for spherical swimmers and being maximal for elongated bodies. Finally, at larger shear rates, the active suspension becomes non-Newtonian Saintillan (2010); Gachelin et al. (2013).

The interactions among swimmers are far from simple, but they can be classified as long- and short-range. The former, of hydrodynamical origin, leads, for instance, to instabilities for pushers Saintillan and Shelley (2008) and to spiral vortices in confinement Wioland et al. (2013). The latter, due to steric and hydrodynamic mechanisms, leads to the formation of polar and nematic phases Cisneros et al. (2007); Sokolov et al. (2007), and flocking transitions Vicsek et al. (1995). Density plays a major role in selecting which interaction dominates. At low concentrations, short-range interactions are negligible, while for moderate and high concentrations they take over the long-range hydrodynamics Cisneros et al. (2007); Drescher et al. (2011). This has been observed as an effective “screening” for very dense suspensions Wensink et al. (2012); Marchetti et al. (2013). Finally, as the long-range hydrodynamic interactions are mediated through the divergence of the stress tensor Saintillan and Shelley (2008), this regime appears whenever the characteristic length of the fluctuations is smaller than the average distance among bacteria. In particular, they can be discarded for homogeneous systems. Finally, long-range hydrodynamic interaction can be included with a mean field description, where the local rheology is governed by the short-range interactions Saintillan and Shelley (2008); Guazzelli and Morris (2011). In this article, we aim to characterize the non-ideal rheology of bacterial suspensions in the semi-dilute regime, hence, considering the effects of short-range interactions.

II Kinetic theory description

To simplify the analysis and render more transparent the rheological response, we consider a bacterial suspension in two spatial dimensions, with conclusions that can directly be extended to three dimensions. Bacteria are described by their position $\mathbf{r}=x\hat{\mathbf{x}}+y\hat{\mathbf{y}}$ and director, which is characterized by a single angle, $\hat{\mathbf{p}}=\cos\theta\hat{\mathbf{x}}+\sin\theta\hat{\mathbf{y}}$. To analyze the bacterial suspension, we use kinetic theory, which was previously used to quantify the viscosity reduction of a dilute suspension of microswimmers Saintillan (2010). In this framework, the object under study is the distribution function $\Psi(\mathbf{r},\theta,t)$, such that $\Psi{\rm d}x{\rm d}y{\rm d}\theta$ gives the number of bacteria with the specified position and director. The suspension is placed in a uniform shear flow, with velocity $\mathbf{u}=u(y)\hat{\mathbf{x}}=\dot{\gamma}y\hat{\mathbf{x}}$, where $\dot{\gamma}$ is the shear rate.

There are three contributions responsible for the change of the distribution function. First, bacteria move as an effect of the imposed flow and their self-propulsion with velocity $V_{0}$, which we assume to be constant and the same for all swimmers,

[TABLE]

and the orientation changes as an effect of the local shear rate, described by the Bretherton–Jeffery equation

[TABLE]

where $\beta$ is the Bretherton coefficient, which for E. coli can be taken as $0.7$, representing its elongated body. These deterministic equations of motion are transcribed to the kinetic equation as advection terms. The second contribution is the rotational diffusion with coefficient $D_{\text{r}}$, which is accounted for with a Fokker–Planck term. Tumbling, which also changes the orientations randomly, gives rise to rotational diffusion at long times, and its effect is absorbed into $D_{\text{r}}$. The explicit inclusion of tumbling with a Lorentz-like term in the kinetic equation Saintillan (2010); Soto (2016) alters only quantitatively the conclusions of this article. Translational diffusion is not considered, because the Brownian diffusivity is much smaller than the Berg diffusivity $D_{\text{Berg}}=V_{0}^{2}/D_{\text{r}}$ Howse et al. (2007), and its inclusion results only in small corrections. Finally, in the semidilute regime, for homogeneous states where long-range hydrodynamic interactions can be neglected, bacteria have binary interactions. These are accounted for with a Boltzmann-like bilinear collision term with both swimmers located at the same position. Considering all these elements, the kinetic equation reads

[TABLE]

where $J$ is the collisional integral. Short-range steric Peruani et al. (2006); Baskaran and Marchetti (2008); Wensink et al. (2012) and hydrodynamic interactions Ishikawa et al. (2006); Llopis and Pagonabarraga (2010); Lushi et al. (2014) align bacteria. For simplicity, we consider that the alignment is total, meaning that after their interaction, both swimmers emerge with the same angle, equal to the average of the incoming ones, resulting in

[TABLE]

as was introduced in Refs. Aranson and Tsimring (2005); Aranson et al. (2007) to describe the alignment of microtubules and bacteria. Here, $g$ quantifies the collision rate. In two dimensions, it scales as $g=\sigma_{2D}V_{0}$, where $\sigma_{2D}$ is the two-dimensional cross section. For E. coli, we can use $V_{0}\approx$20\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{/}\mathrm{s} Berg ([2008](#bib.bib1)) and $\sigma_{2D}\approx$2\text{\,}\mathrm{\SIUnitSymbolMicro m} to estimate $g\approx$40\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\mathrm{/}\mathrm{s}. Also, for *E. coli*, $D_{\text{r}}\approx$0.0125\text{\,}\mathrm{s}^{-1} com .

The total shear stress is the sum of the viscous stress $\Sigma_{xy}^{\text{v}}=\eta\dot{\gamma}$, where $\eta$ is the fluid viscosity, and the particle stress. The latter is the sum of the drag stress for passive suspensions $\Sigma_{xy}^{\text{d}}$ and the active contribution $\Sigma_{xy}^{\text{a}}$ Saintillan (2010). This last contribution, responsible for the effective viscosity reduction, is

[TABLE]

which can be computed once the solution of the kinetic equation is found. Here, $\sigma_{0}$ is the strength of the elementary force dipole, which is negative (positive) for pushers (pullers). For active suspensions, $\Sigma_{xy}^{\text{a}}$ dominates over $\Sigma_{xy}^{\text{d}}$, allowing us to discard the latter.

Equation (3) with $g=0$ (no binary collisions) has been solved in Ref. Saintillan (2010) finding the viscosity reduction described above. The opposite case, with interactions and $\zeta=0$, presents a polar phase above a threshold density Aranson and Tsimring (2005). In both limits, the homogenous case already presents interesting features that could be observed experimentally. In this article, we analyze the homogeneous suspension when both $g$ and $\zeta$ are finite, showing that intriguing rheological features appear.

III Rheological response

Decomposing the distribution function in an angular Fourier series, $\Psi(\theta,t)=\sum_{n}a_{n}(t)e^{in\theta}$, with $a_{-n}=a_{n}^{*}$, the kinetic equation for a homogeneous suspension reads

[TABLE]

where $\text{sinc}(x)=\sin(x)/x$ and we employed the following rescaling to dimensionless variables: $s=D_{\text{r}}t$, $\zeta=\dot{\gamma}/D_{\text{r}}$, and $b_{n}=a_{n}g/D_{\text{r}}$. The global dimensionless density is $\rho=2\pi b_{0}$, which is conserved, implying that $db_{0}/ds=0$. In terms of the Fourier coefficients, the average orientation vector and the swimming contribution to the shear stress are

[TABLE]

The rotational diffusion produces a fast decay of the high $n$ modes. It is possible then to truncate Eq. (6) with few modes, as the rest will take small values. Considering that the swimming contribution to the stress tensor depends on $a_{2}$, we first truncate the series to $n=2$ for the purpose of analytic calculations,

[TABLE]

where we have defined the bifurcation parameter $\epsilon=\rho(4/\pi-1)-1$. When $\zeta=0$, an isotropic to polar transition takes place at $\epsilon=0$, that is, for the critical density $\rho_{c}=\pi/(4-\pi)\approx 3.66$ Aranson and Tsimring (2005). Close to the transition, the growth rate of $b_{1}$ is much smaller than the decay rate of $b_{2}$, and therefore $b_{2}$ is enslaved to $b_{1}$. This analysis leads to the scalings $b_{1}\propto\sqrt{\epsilon}$ and $b_{2}\propto\epsilon$.

The imposed shear flow when $\zeta\neq 0$ induces a macroscopic rotation of $\hat{\mathbf{p}}$ and phase oscillations will appear in $b_{1}$ and $b_{2}$. For $|\epsilon|\ll 1$ a multiple time scale approach is possible, where we define the fast and slow time scales $s_{0}=s$ and $s_{1}=\epsilon s$. Based on the couplings in Eqs. (8) and (9), we propose the following ansatz:

[TABLE]

where $A$, $B$, $F$, $G$, and $H$ depend solely on $s_{1}$. Here, we assume that $\zeta$ is at least one order higher than $\epsilon$ so the two time scales make sense: for the oscillations and for the evolution of the amplitudes.

Using the scalings $A,B\propto\sqrt{\epsilon}$, $F,G\propto\epsilon$, and $H=H_{0}+\epsilon H_{1}$, with $H_{0},H_{1}\propto 1$, Eqs. (8) and (9) are solved order by order in powers of $\epsilon$. To order $\epsilon^{0}$, we obtain

[TABLE]

To order $\epsilon^{1/2}$, it results in the linear system

[TABLE]

from which it is possible to obtain the relation

[TABLE]

and the frequency

[TABLE]

We select the positive root; the negative root is equivalent, with the roles of $A$ and $F$ exchanged with those of $B$ and $G$. To order $\epsilon$,

[TABLE]

Finally, to order $\epsilon^{3/2}$, we obtain the evolution of $A$ in the slow time scale

[TABLE]

in which we have already substituted the previous results.

For $\epsilon<0$ ($\rho<\rho_{c}$), the stationary solution is $A=0$, and the system reaches an isotropic phase where the only contribution to the active stress comes from the imaginary part of $H_{0}+\epsilon H_{1}=H_{0}+\epsilon h_{4}$. On the contrary, for $\epsilon>0$, a polar phase develops with $A\neq 0$. Then, $F$ and $G$ are nonzero, resulting in a shear stress that oscillates in time. For a spherical swimmer ($\beta=0$), $\omega_{0}=\zeta/2$ and $l=0$. Using $\beta=0.7$, which represents E. coli, $\omega_{0}=0.49\zeta+{\textit{O}}(\zeta^{2})$. In this case, $l\approx 0$ and hence $B\approx 0$, meaning that the average orientation vector rotates almost harmonically in the same direction as the vorticity.

To obtain simple analytical expressions, while retaining the main effects, we can neglect $l$. In this case, the stationary solution is $|A|^{2}=\frac{3\epsilon}{8}\frac{1-i\omega_{1}}{h_{1}}$. The latter solution is valid only if it is a real positive number. This condition gives $\omega_{1}=-\text{Im}(h_{1})/\text{Re}(h_{1})$, and therefore the finite value for $|A|^{2}=\frac{3\epsilon}{8\text{Re}(h_{1})}$. In this approximation $G=0$ and from Eq. (7), the active shear stress oscillates as

[TABLE]

with amplitude and mean value

[TABLE]

Having retained $l$ would give finite values for $B$ and $G$, which manifest in nonharmonic terms in the oscillations of $\Sigma_{xy}^{\text{a}}$. Numerical solutions of Eqs. (6), truncating at $b_{\pm 5}$, confirm the qualitative behavior obtained in the multiscale scheme, for small values of $\zeta$. However, a different scenario appears for high shear rates, for which the polar phase disappears continuously, as shown in Fig. 1. This feature appears only when more than two Fourier modes are considered, otherwise the polar phase remains up to arbitrarily large values of $\zeta$. The maximum shear rate for the polar phase depends on the Bretherton parameter, with a phase diagram presented in Fig. 2.

IV Axial alignment

For elongated microorganisms such as Bacillus subtilis or Paramecia, different observations indicate that they interact axially Kemkemer et al. (2000); Kaiser (2003); Ishikawa et al. (2007); Ishikawa and Hota (2006); Zhang et al. (2010); Marchetti et al. (2013), i.e. swimmers with opposite orientations tend to align in an antiparallel way, whereas if they are swimming in a similar direction they align in a parallel manner. As before, to facilitate the calculations we consider the total aligning case. In this situation, the alignment is polar if the pre-collisional relative angle is smaller than $\pi/2$ and apolar if it is greater, which is described by the binary collision operator

[TABLE]

The evolution for each of the Fourier amplitudes is given by

[TABLE]

using the same dimensionless variables as in (6). Note that the evolution of the even modes is decoupled from the odd ones, which is a direct consequence of the nematic symmetry present in the collision rule and the swimmer shape. Equations (23) are solved numerically, truncating at $b_{\pm 10}$.

We define the bifurcation parameter $\nu=\rho(4/\pi-1)-4$, such that for $\nu=0$, there is a transition from the isotropic to nematic phase when $\zeta=0$. This time the critical density is four times bigger than in the polar case, $\rho_{c}^{\text{nematic}}=4\rho_{c}$. When $\zeta\neq 0$, the nematic phase leads to a similar qualitative behavior to the polar case, namely, an oscillating active shear stress. However, important differences appear. First, the average orientation vector is always zero, indicating that no polar order appears. Second, the transition from the ordered to the isotropic phase takes place for very small shear rates. It is necessary to consider large values of $\nu$ (increase the density) to remain in the nematic phase for finite values of $\zeta$, particularly for very slender swimmers (see Fig. 2, bottom). As a consequence that $\nu$ is not necessarily small, there is no full enslaving of the higher modes to the critical one, $b_{\pm 2}$, and the oscillations become non-harmonic, contrary to the polar aligning case, where the oscillations remain roughly harmonic for all parameters. Notably, due to the anharmonicities, the average active shear stress can be positive instead of negative, meaning that pushers interacting axially increase the effective viscosity (see Fig. 1, bottom). Finally, at the transition from the nematic to the isotropic phase, the amplitude of the oscillation vanishes abruptly, and the period diverges with extremely nonharmonic dynamics (see Figs. 1 and 3). These are characteristic features of the limit cycle disappearing via a homoclinic bifurcation Vitt et al. (1996); Strogatz (2018).

V Discussion

In this article we have reviewed the nonideal rheology of bacterial suspensions in the semi-dilute regime. Considering both polar and axial aligning collisions, we have obtained that the short-range interactions lead to an oscillatory active shear stress, in response to a constant shear rate. In both cases we found that there is a competition between the rotation imposed by the shear rate and the short-range interactions. When the shear rate is large, the ordered phase is lost in a continuous or discontinuous transition, for polar or axial interactions, respectively.

To our knowledge, just one experiment has measured with high precision the shear stress in bacterial suspensions López et al. (2015). There, using an E. coli suspension, they observed a large reduction of viscosity (up to a superfluid state) for volume fractions in the range $0.11\%\leq\phi\leq 0.67\%$, which are higher than the critical density predicted in this model for polar aligning interactions ($\phi_{\text{c}}\approx 0.1\%$). These volume fractions correspond to the bifurcation parameter in the range $-0.03\leq\epsilon\leq 4.9$ or $-3.03\leq\nu\leq 1.9$. The experiments used a high-precision Taylor-Couette rheometer subjected to a small shear rate $\dot{\gamma}=$0.04\text{,}\mathrm{s}^{-1} ($\zeta\approx 3.2$). Although these experiments fell well inside the polar phase (see Fig. [2](#S3.F2)), no oscillation in the shear stress was observed. This can be understood by at least three facts. First, the period of the oscillations is $T\sim 2\pi/\dot{\gamma}\approx$150\text{\,}\mathrm{s}, which is larger than the $30\text{\,}\mathrm{s}$ that the experiment lasted and hence there was no time to observe the oscillation. Second, in an extended system, it takes a long time for the different domains to synchronize. Then, the rheometer measures an incoherent sum of oscillatory stresses. Finally, the buildup time of the polar phase scales as $1/(D_{\text{r}}\epsilon)$, ranging between $8\text{\,}\mathrm{s}$ and $100\text{\,}\mathrm{s}$, which limits also the possibility to observe the oscillations in the earlier experiment.

In this article, we have focused on pusher swimmers, for which was predicted a reduction of the viscosity in the dilute limit Saintillan (2010). However, two remarkable effects appear. With axial alignment, the average viscosity increases instead. And the shear stress oscillations involve also an increase of the viscosity during half of the oscillation period. Finally, this implies that pullers should reduce the suspension viscosity under appropriate conditions.

The collision rules employed in this article are totally aligning, and we did not consider any threshold for the relative angle to select which bacteria interact, for example. Nevertheless, we expect the presented results to remain qualitatively the same by relaxing these conditions. For axial aligning collisions, the oscillation period varies quite dramatically near the homoclinic transition, for example, varying the Bretherton coefficient, as shown in Fig. 3, bottom. Then, in a suspension with a diversity of shapes, interesting segregation or synchronization effects can appear. Moreover, a more realistic situation could be accounted for by mixing the two collision rules, with collision factors $g^{\text{polar}}$ and $g^{\text{axial}}$. In this regard, some intuition can be obtained from the study of active polar liquid crystals of slender rods ($\beta\rightarrow 1$) Forest et al. (2014), where these two interactions are dealt with via intermolecular potentials Masao (1981), although there is not a clear relation between the nematic and polar strengths of the potential and the bacterial density. Further studies are needed to understand these effects.

For the three-dimensional case, the distribution function depends on the two spherical angles that parametrize the directors. An expansion of the distribution function in spherical harmonics $Y_{lm}$ would allow one to perform an analogous analysis as in this article. Again, the average orientation vector depends on the $l=1$ modes and the active stress tensor on the $l=2$ modes. Symmetry arguments dictate that the same nonlinear mode couplings will take place, for example, coupling only even modes for the axial alignment case. Hence, although the calculations are much more involved, it is expected that similar rotating polar and nematic phases develop, which will imply that the active stresses will oscillate in time. The eventual end of these phases when increasing the shear rate, in the form of supercritical or homoclinic bifurcations, depends on the explicit expressions for the mode equations and is a subject of future research.

Acknowledgements.

We thank M. Clerc for his comments on the theory of dynamical systems and E. Clement for enlightening discussions on the bacterial rheology. This research was supported by the Fondecyt Grant No. 1180791 and the Millennium Nucleus “Physics of active matter” of the Millennium Scientific Initiative of the Ministry of Economy, Development and Tourism (Chile). M.G. acknowledges the Conicyt PFCHA Magister Nacional Scholarship 2016-22162176.