Counterion-controlled phase equilibria in a charge-regulated polymer solution

We study phase equilibria in a minimal model of charge-regulated polymer solutions. Our model consists of a single polymer species whose charge state arises from protonation-deprotonation processes in the presence of a dissolved acid, whose anions serve as screening counteri-ons. We explicitly account for variability in the polymers’ charge states. Homogeneous equilibria in this model system are characterised by the total concentration of polymers, the concentration of counter-ions and the charge distributions of polymers which can be computed with the help of analytical approximations. We use these analytical results to characterise how parameter values and solution acidity influence equilibrium charge distributions and identify for which regimes uni-modal and multi-modal charge distributions arise. We then study the interplay between charge regulation, solution acidity and phase separation. We find that charge regulation has a significant impact on polymer solubility and allows for non-linear responses to the solution acidity: Re-entrant phase behaviour is possible in response to increasing solution acidity. Moreover, we show that phase separation can yield to the coexistence of local environments characterised by different charge distributions


I. INTRODUCTION
Solutions with charged polymers can demix into polymer-rich phases, also known as condensates.When the condensed phase remains liquid, the process yielding to demixing is known as liquid-liquid phase separation or coacervation.In recent years, the understanding of liquid-liquid phase separation (LLPS) has gained enormous interest because of its putative role in the assembly of macromolecules (mostly proteins and nucleic acids) into membrane-less organelles (also known as biomolecular condensates) in cells [1,2].While polymer physics theories have elucidated several aspects of phase separation in solution, it is not yet fully understood how different molecular mechanisms affect the formation, regulation and properties of biomolecular condensates in cells [2].Challenges relate to the complexity of proteins, that are large heteropolymeric polyelectrolytes, and of the cellular environment which is maintained out of equilibrium and can itself modulate proteins properties and coacervation [2].
Grounded in the seminal work by Flory and Huggins (FH) on phase separation in polymer solutions, the balance between enthalpic and entropic interactions is considered to be the driving force of LLPS.Based on the simplifying assumption of polymers consisting on chemically identical units, Flory and Huggins derived a meanfield model for phase separation in two-components mix-tures.Such a model has proven a useful phenomenological model also to study phase-separation in protein solutions.However, its has limited predictive power, as it misses details on the nature of the intermolecular interactions that contribute to the enthalpic part of the free energy [2,3].
A feature common to proteins is the presence of ionizable groups, that contribute to the electrostatic interactions between proteins [4].Models of polyelectrolyte coacervation are commonly employed to study the role of electrostatic interactions as well as salt in LLPS.The early key paper in the field of polyelectrolyte complexation (also called complex coacervation) remains the work by Voorn and Overbeck from 1957 [5].Extensions of these classical theories that capture the sequencedependence of LLPS driven by proteins with intrinsically disordered domains, as first demonstrated in [6], have employed mean-field theories of polyampholytes as underlying models of proteins.They include the Random Phase Approximation [7,8], as well as Field Theoretic Simulations [9] for the residue specific electrostatic interactions; recent reviews in the modern context are [10][11][12][13].
A limitation of all these approaches is that they assume the charge state on the polymers, such as polyelectrolytes or polyampholytes, to be fixed; in contrast, as shown earlier on by the work of Linderstrøm-Lang [14], the charge state of proteins is in fact regulated by the local environment, such as pH conditions, as well as by interactions between ionizable groups themselves [15,16].
A key process in this context is charge regulation of the polymers or, more generally, chargeable macromolecules in the cellular context [17,18].The charge regulation process is best explained in its most elementary variant which consists in the binding and unbinding of protons, H + , from the water solvent.It is immediately clear that this protonation-deprotonation process goes in hand-inhand with the change of solution pH [19].More involved charge regulation processes are obviously present, e.g. in the binding of dissolved salts in solutions.The effect of charge regulation processes has on phase equilibria has been addressed in several recent papers [20][21][22][23][24][25][26].However, even in simple model systems, the complexity of the interactions yields phase behaviours in multi-parameter spaces which are non-trivial to analyse.This is particularly true due to the highly non-linear free energy terms associated with electrostatic correlation effects, a key feature of liquid-liquid phase separating systems and of fundamental relevance in cell biology [27][28][29][30][31][32][33].
In cell biology, the relation between the phase diagram on the one hand and the charge states of the macromolecules on the other [34] is of particular interest.In this paper, we address this issue on the basis of a 'minimal' model which has essentially two ingredients: a basic formulation of the Voorn-Overbeek theory and the charge regulation mechanism, for which we keep track of the charge state on the polymers following the charge distribution approach developed in [18].Another key novelty that distinguishes our work from previous studies on phase separation and charge regulation processes [20,25] is that we consider the protonation-deprotonation equilibria in solution in the presence of a dissociated acid.The concentration of the counterions due to acid dissociation will turn out to be a key control parameter in our model system.In this way, we are capable to gain insights into the coupling between charge regulation, acidity and phase separation, by linking topological changes in the coexistence curves as well as the related changes in the charge distributions on the polymers.
Our paper is organised as follows.In Section II we introduce our model for the polymer-solvent mixture.Section III covers the results we have obtained from its analysis.Section III A describes its homogeneous equilibrium states, with a focus on how the composition of the mixture affects the polymer charge.Section III B then discusses phase equilibria in our system.Finally, in Section III C we show how phase separation process itself regulates the charge state of the polymers by controlling the local environmental conditions -here acidity.Section IV concludes and provides an outlook to further studies; in particular, we discuss the putative relevance of our results for LLPS in biological systems.Section V contains the Appendices in which the technical results employed in the paper are derived.

II. A MODEL FOR A POLYMER-SOLVENT MIXTURE
Components of the mixture.The building blocks of our model and the charge regulation mechanism it entails are illustrated in Figure 1, respectively.We consider chargeable polymers solvated in water, H 2 O, and a strong acid; here as an example, we consider hydrochloric acid, HCl.Therefore, in solution, we encounter the dissociated ionic species: Cl -and hydronium ions H 3 O + .The polymers are considered as monodisperse with N 1 monomers, of which only a subset of Z monomers carries a protonation site, which can either be positivelycharged (bound state) or neutral (unbound state).We assume that H 3 O + , Cl -and the monomers making up the polymer have the same molecular volume as water, ν, so that the polymers have the molecular volume ν M = N ν ν.As in [18], we assume that polymers with different charge states, z ∈ {0, . . ., Z}, coexist in the mixture; as a result, we have effectively Z + 1 different polymer species in solution.Together with water, chloride and hydronium ions this gives a total of Z + 4 species that we take into account in our mixture.For each species, we denote by φ ω the volume fraction, with ω = (s, +, , z) = (solvent, hydronium ions, chloride ions, charged polymer).The volume fractions must satisfy a no-void condition, which guarantees that at any location space is fully occupied by the mixture: where Furthermore, we assume that our solution is electroneutral so that the net charge density of the mixture has to be zero, The free energy density of a homogeneous mixture.We assume that the mixture is incompressible and kept at a constant temperature T , and describe it by a Helmholtz free energy density f which consists of three contributions, similar to [25], The chemical potentials of the different species in the mixture are then expressed in terms of derivatives of the Helmholtz free energy density f with respect to φ ω ; these conditions are given in detail in Appendix A. The first contribution f 1 in (4) is the standard Flory-Huggins free energy capturing the entropic contributions and an interaction term of water and the solvated polymer For simplicity, we assume the interaction parameter χ to be independent of the charge on the polymers.The where u z is the difference in the internal free energy (nondimensionalised by k B T ) of a polymer with charge z and a neutral one.By neglecting chain connectivity of the polymers, we can see the charged polymer as a mixture of an uncharged polymer and z positive fixed charges.Following [18], we specify u z as In (7), the first contribution represents the energy gain (again non-dimensionalised by k B T ) from occupying an additional site on the polymer by an H + ion.The second term represents an additional contribution from shortrange interactions between occupied binding sites whose strength is controlled by the parameter η.Finally, we have to include the internal entropy to account for the different ways to arrange fixed charges on the binding sites.
The last term f 3 is the Debye-Hückel term, similar to [25], which like our reasoning for f 2 assumes that the charges on the monomers of the polymers can be treated as free ions, where Note that the term κ 2 depends on the sum of all charged molecules multiplied by their valency (as in Eq. ( 6) in [25]).The simplified expression Eq. ( 9) is obtained by applying (3).In Eq. ( 9) the parameter λ = 4π B /a w , where B is the Bjerrum length in water and a w = ν 1/3 is the size of the species in the solution.More realistic models that include polymer connectivity have, e.g., been discussed in [28].However, these include information on the specific location of the charges along the polymer chains.
The charge regulation process.As mentioned in the introduction, models of polymer coacervation commonly assume the charge state on the polymer phase to be fixed.In our framework, this corresponds to assuming that all protonation sites on the polymer are occupied, i.e., imposing in Equations ( 5)-( 9) φ z = 0 for all z = {0, . . ., Z − 1}.We instead assume that charges can reversibly bind to protonation sites according to the reaction where M z represents the polymer with z charges.Then, the charge states of polymers in solution is determined by imposing chemical equilibrium, instead of being prescribed a priori.
Making use of the definition of the Helmholtz free energy (see Section V A) we have that the change in the free energy for each chemical reaction (M z -1 + H 3 O + −− −− M z + H 2 O) occurring in the mixture is given by At chemical equilibrium, Eq. ( 9) must be zero -i.e., the difference in chemical potential of products and reactant of each chemical reaction must be zero.Manipulating Eq. ( 9) we can express µ z in terms of the chemical potential of the counterions, solvent and uncharged polymers: Equation ( 10) can be viewed as an iterative discrete map that, given µ 0 , defines the chemical potential of all charged polymers in terms of µ + and µ s , Using the explicit form of the chemical potential (31) in (11) we arrive at where u z is defined by Equation (7) and Taking the exponential of both sides of (12), we obtain a system of Z+1 linear algebraic equations for the volume fractions φ z ; this can be solved explicitly to obtain an expression for φ z , z = 0, . . ., Z, with where The terms π z indicate the fraction of the total number of polymers in the charged state z as a function of the mixture composition.By definition, their sum must be unity, Z z=0 π z = 1.Inspecting (15), we find that π z can be rewritten in terms of an effective charge regulation free energy, where with The comparison of Eq. (18a) to the definition of u z (see Eq. ( 7)), shows that, in our system, the local composition of the mixture affects the charge regulation process by controlling the energy associated with the protonation/deprotonation of a single binding site.Note the introduction of an effective parameter α eff that includes a composition-dependent correction to the 'bare' linear term in u z .As in [18], we find that the ion concentration in solution, φ + /φ s , affects the effective binding energy.Furthermore, by introducing the Flory-Huggins term in the free-energy, we have that the polymer concentration, φ M , itself affects the binding of ions in solution (see the last term in Eq. (18b)).Using (15) to eliminate φ z (z = 0, . . ., Z) from the definition of free energy density (see ( 4)-( 8)) we obtain the expression for the free energy for an ionic solution with charge regulating polymers where κ = √ 2λφ and we have introduced the variable Q that represents the mean charge of the polymer phase Q = Z z=0 zπ z While we have defined the free energy in terms of the variables φ s , φ + , φ M and φ , the degrees of freedom of the model can be reduced to only two by observing the two constraints (no-void and electro-neutrality) formulated in (1) and (3), that is N Q + φ .These determine φ s and φ + in terms of φ M and φ , albeit, in the case of φ + , only implicitly.

III. RESULTS
In the current work, we focus on the interplay between charge regulation processes and phase separation.Our analysis highlights the key role of parameter η in the equilibrium properties of the system.We, therefore, consider it as a free parameter while fixing the others.
Based on previous works, we set λ = 26.68 [25] and ν ≈ 3.1 × 10 −23 ml [25].The number of monomers in the protein is set to N = 100; of these, we assume that Z = 20 have a H + binding site.We set α = −6.5 so that it is energetically favourable for an individual binding site to be occupied (see Figure 1b).The temperature is fixed to T = 298 K and the Flory parameter to χ = 0.95; the latter value is chosen so that phase separation is observed -even when considering a neutral polymer (see Section III B 2).
A. Analysis of homogeneous equilibrium states.
We first study the properties of homogeneous equilibrium states that arise in our model.We are specifically interested in how the charge distribution of the polymers, π z , depends on the mixture composition, φ M and φ .This is obtained by solving the non-linear system of algebraic equations given by Eqs. ( 1), ( 3) and ( 14)-( 16).Generally, this can not be done analytically and requires numerical approaches.However, we make the following observations.
1.) In the case η = 0 (i.e., independent ion adsorption), the charge distributions π z is binomial, which can be approximated by a Gaussian distribution in z when taking the maximum charge, Z 1; 2.) For z taken as a continuous variable, we can approximate the effective charge regulation free energy u eff as in the limit Z 1 and z ∈ (0, Z) (for the details, see Section V B).In Figure 2, we summarise how the number and location of the local minima of ũeff is controlled by the mixture composition -i.e., the value of the parameter α eff .When ũeff has a single minimum, then we can estimate π z within a saddle-point approximation that we detail in Section V B. We find that for η > −4, we can approximate the charge distribution by a Gaussian distribution whose mean is determined by the minimum of ũeff .
3) The saddle-point approximation is not always valid for η < −4.The breakdown of the saddle-point approximation is due to the appearance of multiple extrema for the function ũeff (see green area in Figure 2) that is reflected in the charge distribution π z having multiple peaks.In this case of failure of the saddle-point approximation, we need to resort to numerical methods of computation.
This general feature of unimodality vs. multimodality of the charge distribution is summarised in Figure 2 which displays the (η, α eff ) diagram.As shown, we can identify two characteristic regimes depending on the value of η: when η > −4, α eff (i.e., the mixture composition), controls the location of the minimizer of ũeff which is always unique; similarly of u eff z .When η < −4, α eff (i.e., the mixture composition), controls both the location and the number of minimizers of ũeff , and likewise of u eff z .We note that transitions between unimodality to multimodality in charge regulating systems had earlier been seen in [18].
We now discuss the three different cases of interest separately in more detail.
By setting η = 0, Eqs. ( 23) are exact and this can be shown without the need of any approximation.Indeed, we have that A can be evaluated explicitly: where Thus the distribution of polymer states, normalised by the total polymer concentration φ M , has the form of a binomial distribution B(Z, p).We can explain the appearance binomial distribution of the charge state of polymers intuitively.When η = 0 there is no correlation of different binding sites; thus the state of each of the Z sites can be treated as an independent Bernoulli random variable with probability of success (i.e., binding) equal to p (see Eq. ( 22)).
When the value of α eff is such that we lie outside the green region in Figure 2, the charge distribution π z is unimodal with most polymers having a charge state similar to z ≈ Q, defined as the unique minimizer of (20).As shown in Section V B, π z can be approximated by a Gaussian whose mean charge Q and standard deviation S , can be written as where p ∈ (0, 1) is implicitly defined by In the case η > −4, S 2 is guaranteed to be positive independently of the value of p ∈ (0, 1).When comparing the exact form of Q and S in the case η = 0 (see Equation ( 21)) and the approximated form for η = 0 (see Equation ( 23)), we find clear parallelisms.When considering η = 0, the model captures the extra energy contributions due to the interaction of the charges on the polymers.Unlike from the case η = 0, this introduces correlation amongst the state of binding sites (occupied or unoccupied) on the same polymer.Nonetheless, we may still interpret p in Equation (23b) as the binding probability for an H + ion to a free binding site.We note that the analogy with the binomial distribution is not exact and difference emerges when comparing the second moments -here the variance S 2 -which explicitly depends on η.When considering states with the same mean charge Q, we have that η > 0 (short-range repulsion) results in a reduction of the variance of the distribution.In contrast, negative values of η yield to wider distributions, i.e., larger values of S. So far, we have considered α eff as a prescribed parameter.However, as illustrated in Equation (23), α eff is determined by the mixture composition -i.e., the values of φ M and φ .The computation of the corresponding concentration diagrams requires solving highly non-linear equations, for which existence and uniqueness of solution may not be guaranteed.Due to the physical constraints in the system (no-void and electroneutrality), homogeneous equilibrium states only exists when φ M and φ satisfy: We can prove that such homogeneous states are unique (see Appendix V C for details).We obtain the solutions numerically via Newton's method and use the approximation to estimate how Q and S vary as a function of the mixture composition.Results for different values of η > −4 are shown in Figure 3.When η is negative (as in Figure 3a), the fully-charged state is the most energetically favourable for the polymers -recall α is also taken to be negative.As a result, whenever the concentration of H + -ions in the mixture exceeds the concentration of the binding sites (i.e., φ > (Z/N )φ M -above the dotted light-blue curve in Figure 3), the polymers will be in a fully-charged stateas Q attains its maximum value (see panel (a)) while S its minimum (see panel (b)).In contrast, when the concentration of H + -ions in the mixture is lower than the concentration of the binding sites (i.e., φ < (Z/N )φ M ), the charges are on average distributed homogeneously between the polymers, Q ≈ N φ /(Zφ M ).This can be shown systematically, by considering the limit α → −∞ when estimating p (results not shown).As η increases it becomes less energetically favourable for H + -ions to bind to the polymers that tend to remain in a less charged state even when φ > (Z/N )φ M .As expected, we find that the largest value of S decreases with η.However, when considering the impact of η on S for a specific mixture composition, there is no general trend.For ion-saturated mixture compositions, S increases with η, while for ion-limiting mixture compositions, S decreases with η.

Multi-modal charge distributions: charge demixing.
We now investigate the equilibrium charge distribution for values of η < −4.As discussed at the beginning of this section, in this regime, the saddle-point approximation breaks down and bimodal charge distributions are expected.
We compute the full charge distribution, {π z }, solving the non-linear algebraic system given by Eqs. ( 1), ( 3) and ( 14)-( 16) using Newton's method with arc-length continuation (used to find good initial guesses for the Newton's step).We conjecture that {π z } is still uniquely defined even when we are in regimes for which the charge distri- bution has multiple peaks (i.e., when we enter the green area in Figure 2); this is strongly supported by our numerical investigation but an analytical proof of the result is beyond the scope of this work.
The results are shown in Figure 4 in which we compare the homogeneous equilibrium states for η = −7 (left column), η = −5 (middle column) and η = −2 (right column).Interestingly, we find that Q is almost insensitive to changes in η (recall that here α = −6.5 0); both below and above the H + -saturation curve the mean charge is not affected by increasing of the short-range attractions between bounded charges (i.e., moving from right to left in Figure 4).In contrast, the behaviour of the standard deviation S changes significantly with η; particularly for mixture compositions below the saturation curve.Overall, we find that the more negative η, the larger the maximum value of S. When η < −4 (see first and middle column in Figure 4), large values of the variance S are attained by allowing charges to be distributed unevenly between polymers -i.e., π z has a bimodal profile (see panels (g) and (h) in Figure 4).For values of η near the critical threshold η = −5 (see panel (h)), we find broad distributions, with polymers in all charge states present in the solution.In this case, the peaks in the distributions occur away from Q (see vertical red line) suggesting that most polymers have a charge state that deviates from the mean.As we take η −4 (see panel (g)), π z becomes more skewed towards the ex-treme states, z = 0 and z = 1, and the large values of S are due to the differential partitioning of the charges rather than π z having a broader support.This is because the intermediate charge states, z ≈ Z/2, become energetically unfavourable and most polymers exist either in a poorly-charged (z ≈ 0) or in a highly-charged state (z ≈ 1).In this regime, changes in the mixture composition only impact the relative fraction of the polymers in poorly-charged and in highly-charged states thus allowing Q to attain all values in the interval [0, Z].From this point of view, the model could be approximated by a two-population model: either neutral or fully charged polymers that coexist under proper conditions.This is similar to the approach adopted in [19].

B. Demixing in solutions of charged polymers
In the previous section we have discussed how charge regulation affects the homogeneous equilibrium states of the mixture.In particular, we find that the mixture composition modulates the equilibrium charge distribution.Due to the physical constraints on the volume fractions -i.e., no-void and electro-neutrality -at equilibrium the mixture composition is well-defined by the volume fraction of two species.Here we have chosen: the total volume fractions of polymers, φ M and counterions, φ .
We now investigate how charge regulation impacts the solubility of charged polymers.The calculation of the phase diagrams follows standard procedures -details are given in Appendix V D. We denote by φ I ω and φ II ω the volume fraction of species in the dilute (i.e., polymer depleted) and condensed (i.e., polymer rich) phases, respectively.Importantly, in constructing the phase diagrams we allow the ions to be distributed asymmetrically between the dilute and condensed phases.As a result, the tie-lines (i.e., the curve connecting coexisting states) can have non-zero gradients.This leads to the mean electrostatic potential being different in the dilute (ψ I ) and condensed (ψ II ) phases.The difference ∆ψ = ψ II − ψ I is known as the Galvani potential [32].For any value of the model parameters the phase diagrams are practically computed in Julia using the BifurcationKit package [35] for numerical continuation.
As mentioned in Section II, most models of phase separation assume that the charges on the polymers are fixed.In order to highlight the role of charge regulation in phase separation, we first investigate demixing for a solution of polymers with a fixed charge, Z.While in the charge regulation (CR) model the charge distribution, π z is obtained by minizing the free energy f (see Eqs (4)-( 9)), in a fixed charge (FC) model, π z is prescribed via a delta function 4)-( 9)), we obtain the free energy for the FC model, f FC , as 25) where κ = √ 2λφ (as before) and u Z is defined as in (7).As before, the system must also satisfy the electro-neutral and no-void constraints (see Eqs. ( 1) and ( 3)).
1. Phase diagrams for macromolecules with a fixed charge.
In Figure 5, we present the phase diagram for increasing values of the fixed charge on the macromolecules, Z.In these diagrams, regions of mixing and demixing are separated by the binodal (or coexistence) curves.Along the binodal, we highlight the gradient of the tie-lines: positive gradients (in red) indicate the counterions concentration is higher in the condensed phase (II); in contrast, negative gradients (in blue) imply counterions accumulate in the dilute phase (I).We note that, besides the constraints (24), in the fixed charge model, the electroneutrality condition also requires φ > Z/N φ M .
Starting from the case of neutral polymers (see Figure 5a), we recover a coexistence curve analogous to the one obtained in previous works on coacervates [11,32].
Here the region of demixing is enclosed by a single open curve (the bimodal) and a unique critical point (highlighted in red) exists.Furthermore, the tie-lines have a negative gradient, suggesting that more counterions accumulate in the dilute instead of the condensed phase.The gradient steepens near the critical point, while tielines are almost horizontal when the counterions are dilute (φ 1).As we increase Z the fixed charge on the polymers (see Figure 5b-5c), the demixing region is affected only for small values of φ ; this is primarily due to intersection of the bimodal curve with the boundary of the feasibility region (φ = Z/N φ M ).Since the latter curve has a positive gradient, this enforces the tie-lines to change their orientation as they approach the boundary of the feasibility region.If we increase the charge on the polymers even further (see Figure 5d-Figure 5e), we find the demixing region shrinks and its topology changes into a closed-loop with the emergence of a second critical point.We find also a complete inversion in the slope of the tie-lines compared to the neutral case.If we were to increase Z even further, the miscibility gap will disappear (results not shown).
We conclude that overall fixed charges reduce the solubility of polymers in solution.
In Figure 6, we illustrate the characteristic topologies of the phase diagram for charge-regulating polymers for different values of η.In these diagrams, regions of mixing and demixing are separated by the binodal curves.In Figures 6a-6d, we depict along the bimodal the mean charge on the polymers, Q.In Figures 6e-6h, we illustrate the same phase diagrams but highlight along the binodal the gradient of the tie-lines (see grey curves).Interestingly, we find that the phase diagrams can be significantly different from each other depending on the value of the charge regulation parameter η.In particular, we find that, for strong short-range repulsion between occupied binding sites, i.e., α+η large and negative (first and second column in Figure 6), the phase diagram presents two disconnected regions of demixing -namely A and B in Figure 6a)-which are enclosed in the demixing region obtained for neutral polymers (see shaded area in Figure 6a).The demixing region A in Figure 6a lies above the H + saturation curve (see Figure 3a and related discussion) and the polymers effectively behave as having a fixed charge of Z = 20.When comparing region A in Figure 6a (or Figure 6b) and the demixing region in Figure 5e, the two overlap exactly.In contrast, the demixing region B lies fully or partially below the saturation curve.The boundary of this region is delimited by coexisting phases that differ both in the local amount of polymers as well as in their charge state -as highlighted by the variation in the mean charge Q.The implication of these results will be investigated in Section III C. As the value of η increases (i.e., it is less favourable for polymers to be in a fully charged state), the two disconnected regions merge and a single demixing region persists (see Figure 6c).Eventually, for η sufficiently positive, the phase diagram converges to the one of neutral polymers (see Figure 6d).
Interestingly, when comparing phase diagrams with two demixing regions, we find that the tie lines have always a positive gradient -i.e., the concentration of counterions is lower in the dilute (I) instead of condensed phase (II).In contrast, for the phase diagrams with a single demixing region, we observe different trends in the tie-lines: (e)-(g) always a positive gradient; (h) a mix of tie-lines with positive and negative gradients in the proximity of the critical point.
Overall, we find that, similarly to fixed charges, the presence of charge-regulating binding sites lowers the demixing tendency of polymers (compared to the neutral case -see shaded area in Figure 6).Nonetheless, we find that charge regulation mechanisms, unlike fixed charges, yield more complex topologies of the phase diagrams.As investigated in the next section, this gives rise to nonlinear dependencies between the polymer solubility as a function of the solution acidity.

The impact of counterions on polymer solubility.
Recent studies have focused on studying how chemical properties of salt ions (such as the counterion radii) affect the solubility of charged polymers with fixed charges [36].Their theoretical results, for a system of polyelectrolytes in a solvent with salt (i.e.positive and negative mobile ions), show non-monotonic salt concentration dependence where salting-out at low salt concentrations is due to ionic screening.In the high salt concentration regime, the macromolecules remain in the salting-out regime for small ions but change to a salting-in regime for larger ions.They conclude that the solubility at high salt concentrations is determined by the competition between the solvation energy and the (translational) entropy of ions, addressing the intensely discussed problem of salt effects in LLPS of protein solutions, such as re-entrant phase transitions shown experimentally in [37,38].
Here, we are interested analogously in studying the impact of counterions (or solution acidity) on the solubility of charged polymers.We define the solubility, ω = ω(φ ), of a charged polymer for a given counterion concentration φ , as the minimum value of the equilibrium volume fraction on the binodal curves (see schematic drawing in Figure 7a).Our definition is analogous to the one used in [36], but corrected for the fact that, in our model, multiple coexistence curves may exist.
As shown in Figure 7b, we find that for neutral molecules the solubility increases with counterion concentration (see purple curve).In contrast, when considering polymers with fixed charge, ω, has a non-monotonic profile which agrees with the results obtained in [36] (for relatively large salt ions), despite our simpler approximation of electrostatic fluctuations.At low counterion concentrations, the solubility of the polymers decreases with φ .This trend -which we referred to as conterion-out behaviour -is considered to be universal for all ions at low ionic concentrations and is explained by the fact that the counterions are able to screen the charge on the poly- mers and hence reduce the Coulomb repulsion between the polymers.In contrast, at higher counterion concentrations, the solubility increases with φ -the conterionin effect.This can be explained by the dominant contribution of the entropy of mixing associated with the ions over charge-screening effects, which favours the miscibility of the solution, very similar to the properties of the system studied in [36].
As shown in Figure 7c, solubility curves of chargeregulating polymers more complex trends.When η ≤ 0, we find that the solubility curve can be split into three regimes: acid-in at extremely low counterion concentrations; acid-out for intermediate-to-low counterion concentrations; counterion-in at high counterion concentrations.Note that in the transition between the counterion-in at extremely low φ to counterion-out behaviour for low φ , the solubility curve is not smooth.Jumps in ω and ω is a signature of the presence and merging of the two disconnected demixing regions (see the curves with η ≤ 0 in Figure 7c).For larger values of η (see red curve in Figure 7c), corresponding to the scenario where binding of the ions to the monomers is unfavourable, we recover a monotonic solubility curve as for neutral macro-molecules: consistent counterion-in behaviour (independently of φ ).
The transition in the sign of the first derivative from ω > 0 to ω < 0 is a signature of another important feature of the phase diagrams in Figure 6a-6c: counteriondriven re-entrant phase separation.Specifically, when short-range repulsion are not too strong (see Figures 6a-6c), the system exhibits re-entrant behaviour when varying the concentration of counterions, φ .In other words, there are values of φ M that lie in the demixing region at very low and high values of φ but not for intermediate (or very high) concentrations of counterions.
C. Regulation of the charge distribution via phase separation.
In the previous section, we have shown how charge regulation affects phase separation in solutions of charged polymers.Conversely, in this section, we are interested in how phase separation itself regulates polymer charge in solution.In order to investigate this aspect, we consider a standard quenching experiment where we drive the system to phase separate by controlling the acidity of the solution (i.e., decreasing φ ).Specifically, we start from a homogeneous mixture (O) with composition φ O M = 0.2 and φ O = 0.04; this is then perturbed by decreasing the acid volume fraction to φ = 0.023.When considering spatially homogeneous equilibria, at any location in space the charge distribution of the polymer phase is the same.However, this is not guaranteed when considering a demixed solution consisting of a dilute (I) and condensed (II) phase.In this case, we denote by π I z and π II z the charge of polymers in each of the two phases.When considering the solution as a whole, the charge distribution on the polymers can be expressed as the weighted average of π I z and π II z : where γ is the fraction of the total volume of the solution occupied by the dilute phase (I) in the quenched state (O ).The value of γ is constrained by the conservation of the total concentration of any of the species in the solution; without loss of generality we here consider the conservation of the polymer molecules to obtain: As a reference case, we test the protocol on a solution of non-hydrophobic polymers that do not phase separate -see Figure 8A.In this case, decreasing the acid concentration in the solution (i.e., equivalent to decreasing φ ) does not lead to phase separation.Yet, it significantly affects the polymer charge distribution (compare Figures 8b and 8c), leading to discharging of the polymer binding sites.In Figure 8B, we consider the same ideal protocol applied to a solution with hydrophobic polymers that phase separates in solution when decreasing the acid volume fraction (see Figure 8d).As shown in Figure 8e, the initial charge distribution on the polymers is similar to the one observed on non-hydrophobic polymers (compare with Figure 8b).Upon quenching, the solution phase separates -state (O ) in Figure 8B.Polymers in the dilute phase remain highly charged (see Figure 8g) as in the initial state (O), whereas polymers in the condensed phase partially discharge (see Figure 8f) as in the case of non-hydrophobic polymers (see Figure 8c).When considering the overall solution, the different charge distribution in the two phases is reflected in the charge distribution π O having multiple peaks -see Figure 8h.By controlling the mixture properties locally -here solution acidity -phase separation creates two environments: the condensed phase where the charge distribution on polymers is highly sensitive to changes in the solution acidity, and the dilute phase where the charge distribution is robust to the changes in acidity.As a result, phase separation allows spatial confinement of polymers with a specific charge state.Note that, after quenching, in Figure 8A, polymers with intermediate charge appear homogeneously in the solution, while these are only localised in the condensed phase in Figure 8B.The possible functional implication of these findings in the context of biomolecular condensates will be discussed in the next section.

IV. CONCLUSIONS AND DISCUSSION
In this work, we considered a minimal model to investigate the interplay of phase separation and charge regulation.For this, we introduced in Section II a system of chargeable polymers, whose charge state is regulated by protonation/ deprotonation processes, in a water-acid solution.
In Section III A, we established the homogeneous equilibria states of the system focusing on how the mixture composition -i.e., the concentration of the polymers, φ M , and the counterions, φ -affects the polymer charge distribution.In doing so, we employed analytical findings which highlighted the key role of the parameter η, describing bounded charge interactions, in determining the properties of equilibrium charge distribution.
Our key findings are: For η = 0, the charge distributions in homogeneous states of the system simplified considerably and it can be found to follow a binomial distribution.For η = 0 we showed that by approximating the charge interaction as a continuous function, we can approximate the charge distribution by a Gaussian distribution for a continuous variable in the limit of a large number of charges, as we derive within a saddle-point approximation.Our analysis yielded that for η < −4, this approximation ceases to be generally valid, as multimodal distributions can arise depending on the mixture composition.
In Section III B, we unfolded how charge regulation processes affect phase diagrams of polymer solutions.To do so, we first characterised phase diagrams assuming a fixed charge on the polymers; we then investigated how the topology changes by introducing charge regulation mechanisms.We found that charge regulation processes can affect the phase diagram topology in a nontrivial manner: upon decreasing η, we observed that the usual demixing region undergoes a change of its topology, in which a closed-loop region branches off from the original demixing region -which persists.This contrasts with the phase diagrams of polymers with fixed charges, where at most one demixing region exists.The complex topology of the phase diagram is reflected in the relation between counterions concentration and polymer solubility (in short solubility curves) in an acid-water solution.We find that charge regulation mechanisms have a prominent signature: depending on the charge-interaction parameter η, due to the re-entrant phase behaviour induced by charge regulation, the solubility curve can exhibit a pronounced jump.This might be a relevant experimental signature for the charge-regulation induced transition in the topology of coexistence curves.
In addition, our results show that charge regulation has an important impact on the partitioning of counterions and thus the gradients of the corresponding tie-lines, which is further complicated by the existence of multimodal equilibrium states.These findings add to the discussion on salt-partitioning in the complex coacervation of polyelectrolyte [12].Different theoretical frameworks -e.g., random phase approximation (RPA) which includes connectivity of the polyelectrolyte, and Liquid-State theories-have been proposed to explain salt partitioning and tie-lines gradients.Identifying the physical reasons for salt partitioning amongst the increasing number of candidate theories remains an open problem and an active area of research [12].Our results suggest that allowing for charge regulation can also affect tie-lines gradient thus adding yet another layer to this discussion.
In the last section, Section III C, we investigated the effect phase separation has on the charge distribution.By discussing an experimental scenario in which the concentration of counterions in the polymer solution is changed, we demonstrated that phase separation can create local environments with very different charge distributions in the dilute and condensed phases, which in addition are either very similar or quite different from the initial state.Our findings highlight how charge regulation mechanisms can have a significant role in the response of polymer solutions to changes in the physical environment, by introducing a complex coupling between processes occurring at the micro-scale (protonation/deprotonation) and meso-scale (phase separation).Interestingly, similar nonlinear effects -like e.g.re-entrant phase behaviour-have been recently discussed by Jacobs et al. in a seemingly unrelated system, in which a different molecular process -polymer self-assembly -is discussed in conjunction with phase separation [39].This suggests re-entrant behaviour might be a general feature of systems where phase separation is coupled to a molecular mechanism (such as charge regulation or self-assembly).
On the one hand, charge regulation controls the sensitivity of polymer solutions to environmental changes by allowing for so-called re-entrant demixing behaviour.We expect this non-linear dependence of phase separation on environmental cues to be fundamental in a range of applications to soft matter science, such as in the design of responsive materials, as well as in LLPS of proteins.Salt-induced re-entrant phase separation has been observed for proteins that undergo LLPS in the high-salt regime [37], including the intensively investigated protein FUS.These observations, together with our and previous theoretical works [19] show the impact of the environment as a driving force of LLPS and adds a further important mechanism to the widely discussed sequencedependence LLPS of intrinsically disordered proteins.
Conversely, we also find that by affecting the local environment the polymers are in, phase separation itself can regulate the polymer charge state by allowing to spatially confine polymers in a specific charged state -hence increasing their local volume fraction.This can have important consequences when considering polymers interacting with additional chemical agents, whereby their interactions may be mediated by the polymer charge state.This is the case in the cell cytoplasm.From this point of view, phase separation in cells might function as a regulator of cellular responses to the environment by controlling both the location, as well as the charge state of proteins.
Despite its simplicity, our model yields a rich and interesting range of behaviours that hint at the importance of charge regulation mechanisms in the formation and properties of condensates.There is therefore scope to extend our theory to investigate whether our findings have relevance to LLPS in cells.This requires extending our model to account for the complexity of biological macromolecules -such as RNA and proteins.For example, in this work, we have assumed the interaction parameter χ to be independent of the polymer charge state.However, when considering short-range interactions e.g. between polymer chains, these are known to be chargedependent.Therefore, a further natural extension of this work would be to analyse the scenario in which χ is considered a function of the charge state z.This would result in the mixture composition influencing not only the association-dissociation energy parameter, α, but also higher-order interactions between binding sites.
Overall, our results reveal that, even in the simplest system consisting of one polymer species whose charge states undergo a protonation/deprotonation process, the interplay between phase separation and charge regulation mechanisms governs the response of polymer mixtures to environmental changes.phase of this project.GLC is supported by the UK Engineering and Physical Sciences Research Council (EP-SRC), grant number EP/W524335/1.

V. APPENDIX A. Derivation of the chemical potential condition
We consider an incompressible mixture in the (T , V , N )-ensemble with temperature T , volume V and particle numbers N ω -where ω ∈ Ω.At equilibrium, such a system minimises the Helmholtz free energy, F = F (T, V, {N ω } ω∈Ω ).From Euler's relation, it follows that where p is the pressure and µ ω are the chemical potentials of the different components of the mixture.Incompressibility of the mixture implies that the molecular volume ν ω of each component of the mixture is constant; as a result, the volume of the mixture can not be taken as an independent variable but rather as a function of the particles numbers: V = ω∈Ω N ω ν ω .Differentiation of F with respect to particle number leads to the chemical potential condition Transforming now to the Helmholtz free energy density f (T, {φ ω } ω∈Ω ) = F/V with the volume fractions φ ω = N ω ν ω /V we obtain, performing the necessary differentiations, (30) Applying this general relation to our mixture we find that the chemical potential of the free ions (µ + , µ ), solvent (µ s ) and z-charged polymers (µ z ) are given by The expression for Σ arises from those terms in the round brackets in (30) that do not cancel out which is only the case for contributions from f 2 and f 3 .Making use of the no-void condition for f 2 and the electroneutrality condition for f 3 one finds

B. Saddle-point approximation
When considering η = 0, the charge distribution, {π z } Z z=0 , is binomial.It is well known that a general binomial distribution, B(Z, p), is well-approximated by a Gaussian distribution with the same mean and standard deviation, in the limit Z 1 -provided p is bounded away from its extreme values 0 and 1.Here we show that a saddle-point approximation to the charge distribution is possible provided that η > −4, guaranteeing that the charge distribution has a unique maximum.
Substituting the definition of u z (see (7)) into ( 14), we obtain that π z reads where the α eff is as defined in (18b).In what follows, we want to approximate the distribution (36) by a Gaussian distribution centred at its mean value Q = Z z=0 zπ z under the assumption that Z 1.For our approximation to be valid, the mean charge needs to be sufficiently far from its extreme values, i.e., Q 0 and Z − Q 0. Since Z 1, we rewrite the discrete charge distribution (36), as a continuous probability distribution for the continuous variable z ∈ [0, Z].First, we approximate the binomial coefficient by using Stirling's series: Using (37), we find Following the continuous approximation, we can write (36) but considering z ∈ (0, Z) as a continuous distribution: where By computing the second derivative of ũeff , it is apparent that for η > −4 the function ũeff is convex for ω ∈ (0, 1).This guarantees that there exists a unique minimum, p ∈ (0, 1).As discussed in Section III A, we can interpret p as an effective binding probability of ions to the polymer that we have defined in the main text.An implicit definition for p can be obtained by solving ũ eff (p) = 0: When Z 1, the mass of the normalisation integral (see first factor in Equation (39a)) will be localised around the stationary point, p, and standard techniques, such as Laplace's method can be applied: Substituting the above into Equation (39a) and expanding around the stationary point (Zp), we obtain a Gaussian distribution: FIG. 9: Saddle-point approximation.Plots comparing the exact discrete distribution (36) (histogram) and its approximation obtained via the saddle-point approach (43) (red curve).Different panels corresponds to different choices of the parameter α eff .
In Figure 9, we compare the approximated distribution (43) with the real distribution (36) for different values of α eff .We find good agreement between the two.Nonetheless, discrepancies emerge when considering |α eff | 1 when the maximum of the distribution shifts towards the boundary of the domain: for α eff large and negative, the maximum ≈ Z while for α eff large and positive ≈ 0. This discrepancy is to be expected since for the approximation to hold we must assume p is bounded away the extreme values 0 and 1.
C. Unimodal distribution: domain of physicality.
In this section, we outline results on the existence and uniqueness of the effective binding probability p (see (23b)) assuming η > −4 for any values of φ M ∈ (0, 1) and φ ∈ (0, 1) which are physically allowed.
By substituting ( 23) into ( 1)-( 3), we find that p is implicitly defined by the non-linear algebraic equation Π η (p) = 0, where The form of Π η is obtained starting from ( 40) and (18b) by first eliminating φ + via (3) in the form φ + = φ − Zp/N φ M , and by finally using (1) in the form φ s = 1 + (Zp/N − 1)φ M − 2φ to eliminate the remaining dependence on φ s .Note that when setting η = 0, Π 0 reduces to a quadratic equation for p 0 that can be solved explicitly: Nonetheless, we consider the more general case η > −4, and prove that there exists at most one root p for the function Π η in the interval (0, 1]; conditions for existence are then discussed.We here exclude 0 since p = 0 refers to the critical case where no counter-ions are present in the solution: Π η (0) = −e −α−χφ M φ = 0 only if φ = 0. First, we note that, for x ∈ (0, 1], the first term on the right hand side in (44) is always non-negative (since from the no-void condition φ s + φ + = 1 − φ − φ M ≥ 0).The sign of the second term instead depends on the value of (Zx/N )φ M − φ .Given that p is a root of Π η , if that exists, we must have that Zp/N φ M − φ < 0 (which guarantees that φ + > 0).
We now consider the value of the first derivative of Π η and evaluate it at one of its possible roots: where we have used the fact that Π η (p) = 0.It is apparent that the first term in (46) is positive and so is the second term, since, as discussed above, we must have that φ − Z/N pφ M > 0. This implies that for any root of Π η , p ∈ (0, 1], the derivative Π η (p) > 0.
a. Uniqueness.Since the function Π η is analytic and Π η (p) > 0, we conclude that if p exists this must be unique.Otherwise, there would exists a root p ∈ (0, 1) such that Π η (p) ≤ 0.
b. Existence.Generally, the existence of p is not guaranteed.Since Π η (0) < 0, the only conditions for the existence of p is that Π η (1) ≥ 0: c. Domain of physicality.Summarising the results above, we find two inequality constraints that homogeneous equilibria exist (i.e, p is well-defined) provided that: 1 − φ M − φ > 0, (48a) However, for the equilibria to be physical meaning, we must have that the corresponding volume fractions φ + and φ s are positive and less than one.Conditions (48) are sufficient to guarantee this is the case.

D. Two-phase coexistence conditions
In this section, we derive the coexistence conditions used to compute the phase diagrams presented in Section III B. We start by considering an initially homogeneous mixture of the Z+4 species that has been quenched into the unstable regime, just before separates into two phases.Each of the emerging phases are homogeneous with a unique composition, characterised by the composition vectors φ I ω and φ II ω .In the demixed state, the conditions for the coexistence of two phases are µ({φ I ω } ω∈Ω ) = µ({φ II ω } ω∈Ω ). (49) These are Z + 4 conditions for 2(Z + 4) variables, leaving Z + 4 degrees of freedom.For the charge regulation (CR) model, we assume each phase is in chemical equilibrium, which imposes the chemical potentials in each of the two phases to satisfy Eq. ( 10), or Z restrictions each.When considering the fixed charge (FC) model, the system is constrained by imposing the charge distribution π z = δ(z − Z) in both phases; also in the latter case, this leads to 2Z restrictions.However, due to the equality of chemical potentials between phases, we only need to impose these Z conditions on one phase (for the other they are then implied).So we have 4 degrees of freedom left.We also have to satisfy electroneutrality and no-void in each phase, which removes all four remaining degrees of freedom.As a result, we lack one degree of freedom required to match those of the initial homogeneous mixture prior to demixing.This problem is frequently addressed by adding an additional contribution to the chemical potential µ ω in (31) for each of the charged species, giving rise to the electrochemical potential, μω = µ ω + z ω eψ, (50 where ψ is the Galvani potential.These electrochemical potentials are then equated instead of the chemical potentials.In a homogeneous system, the Galvani potential ψ is constant and hence can be eliminated by setting it to zero, but in a non-homogeneous e.g.demixed system it is usually not.We then have two different values for ψ and the difference between the two remains as the previously missing additional degree of freedom. Here, we proceed differently to motivate the introduction of a Galvani potential and describe the phase separation as a minimisation problem.We again consider a system with two coexisting phases and a total volume V = 1 (without loss of generality), split into two subsystems I and II of volume γ and 1 − γ, respectively, with 0 < γ < 1.Each subsystem is occupied by a single, in itself homogeneous phase described by the variables φ I = φ I ω ω∈Ω and φ II = φ II ω ω∈Ω , respectively.The total free energy of the demixed system is then given by (51) In a system without chemical reactions, each species is individually subject to mass conservation and we would minimise F under these Z + 4 constraints to find the equilibrium of the system.With chemical reactions, a smaller number of quantities are conserved, and these quantities need to be determined in an additional step prior to the formulation of the minimisation problem.For this purpose, note that the total number of molecules of species ω present is given by For to be conserved, the vector a = (a ω ) ω∈Ω has to satisfy where S is the stoichiometric matrix (with Z + 4 rows and Z columns), that is, the rows of its transpose are the stoichiometric coefficients of the chemical reactions.To write out this matrix, we assume that the indices ω are ordered as z = 0, 1, . .
These are the equations we solve using bifurcation packages as described in the main text; f * is substituted by either f CR (see (19)) or f FC (see (25)) depending on the formulation of the model of interest.

FIG. 2 :
FIG. 2: Composition-dependent charged states.Parameter diagram for the charge distribution of homogeneous states as a function of α eff and η, obtained by computing the extrema of ũeff (see (20)).The insets show ũeff for specific values of α eff and η.In the green region ũeff has two minima; outside this region a unique minimum exists and its position is indicated by the colorbar above the diagram.The change of the effective binding energy parameter α eff (red path in the (φ M ,φ )-plane on the left corresponds to moving along a horizontal line in the (α eff ,η)plane).

FIG. 3 :
FIG. 3: Composition dependence of mean charge and standard deviation.(a)-(d) Series of surface plots illustrating how the mean charge, Q, depends on the local composition of the mixture (φ M , φ ), for different values of the parameter η -from left to right: η = −2; η = 0; η = 2 and η = 5 (short-range repulsion between bounded charges).(e)-(g) Same as panels (a)-(d) but illustrating the computed standard deviation, S .The dotted light blue lines indicate the salt concentration at which the concentration of H + ions in solution equilibrates the concentration of binding sites, i.e. φ = Zφ M /N .Other parameters are set to default values given at the start of Section III.

FIG. 4 :
FIG. 4: Composition dependence of the mean charge and standard deviation for η < 0. Series of surface plots illustrating how (a)-(c) Q and (d)-(f), depends on the local composition of the mixture (φ M , φ ), for different values of the parameter η: (left column) η = −7; (middle column) η = −5; and (right column) η = −2 (same as Figure 3a).The dotted light blue lines indicate the salt concentration at which the concentration of H + ions in solution equilibrates the concentration of binding sites, i.e. φ = Zφ M /N .(g)-(i) Plots of the charge distribution, π z (see (14)) for specific values of (φ M , φ ) (see white dots in panels (a)-(f)); the red vertical lines indicate the mean of the distribution, Q.Other parameters are set to default values given at the start of Section III.

FIG. 5 :FIG. 6 :
FIG. 5: Phase diagram topologies for polymers with fixed charges.In the different panels, the following fixed charge values Z have been chosen: (a) Z = 0, (b) Z = 5, (c) Z = 10, (d) Z = 15 and (e) Z = 20.The colour scale indicates the gradient of the tie-lines while tie-lines connecting coexisting states are indicated in light grey.The area of the (φ M ,φ ) space that are unphysical for our model (i.e., electroneutrality is not satisfied) are shadowed in grey.Critical points at which the two coexisting phases become indistinguishable are denoted by the red circles.Other parameters are set to default values given at the start of Section III.

FIG. 7 :
FIG. 7: Counterion effect on polymer solubility.(a) Schematic showing how the solubility, ω, is computed starting from the phase diagrams in Section III B (details in the main text).(b) Solubility ω as a function of the counterion concentration for polymers with different fixed charges (same parameters as in Figure 4a).(c) Solubility ω as a function of the counterion concentrations for charge regulating polymers (same parameters as in Figures 6).

FIG. 8 :
FIG. 8: Phase separation as a charge regulation mechanism.Effect of quenching the solution by decreasing the concentration of counterions.We consider two cases: (A) a solution with non-hydrophobic polymers (χ = 0); (B) a solution of hydrophobic polymers that phase separates upon quenching (χ = 0.95).(a) Phase diagram for the case χ = 0 (no demixing).(b) Charge distribution in the initial mixed state (O).(c) Charge distribution, π z , after the quenching -i.e., homogeneous mixed state (O ).(d) Phase diagram for the case χ = 0.95 (same as in Figure 6c).Decreasing the concentration of counterions drives demixing of the solution into a dilute (I) and condensed (II) phase which are determined by the tie lines.The final state (O ) is a demixed solution where 67% of the volume is occupied by the dilute phase while 33% by the condensed phase.(e) Charge distribution in the initial mixed state (O).(f) Local charge distribution for polymer in the condensed phase (II).(g) Local charge distribution for polymer in the condensed phase (I).(h) Overall charge distribution for the demixed mixture (O ).Parameter values are set to default values and η = −2 (as in Figure 4b).