Sheeba Jem

Richard H Henchman
email: henchman@manchester.ac.uk

Prediction and interpretation of the hydration entropies of monovalent cations and anions

Keywords: solutions, free energy, statistical mechanics, water, ion

The hydration entropy of the monovalent cations and anions are calculated directly from molecular dynamics simulations of the hydrated ion and bulk water using theory previously applied to the hydration of noble gases [Irudayam and Henchman, J. Phys.: Cond. Matter, 2010, 22, 284108]. Extensions are included to account for differential hydrogen-bonding of first-shell waters with themselves, the ion, and bulk water. The entropies, enthalpies and Gibbs energies agree reasonably with simulation and experiment when the effect of force field is taken into account. The anions' entropy losses are mostly vibrational and librational, consistent with their stronger hydration. The cations' entropy losses are mostly orientational which implies fewer hydrogen-bond arrangements because the cation substantially inhibits the ability of surrounding water molecules to accept hydrogen bonds. Owing to the many entropy terms and different decompositions, it is shown that the terms of kosmotropes and chaotropes must be appropriately applied so as not to lead to contradictions. It is also proposed that the number of hydrogen bond arrangements help explain the ordering in Hofmeister series of ions whereby anions increase this number but cations decrease it.

Introduction

Much knowledge has been accrued about the properties of ions in aqueous solutions. [START_REF] Franks | Water: a comprehensive treatise[END_REF][START_REF] Ohtaki | [END_REF][3][4][START_REF] Kunz | Specific Ion Effects[END_REF][START_REF] Zhang | [END_REF] However, one quantity for which our understanding is incomplete is the solution's entropy. We might be able to measure or calculate entropy changes for processes such as the hydration of ions but deeper insights come from knowing the entropy itself and how it is distributed between all the constituent molecules. This makes imperative the development of new methods with these capabilities. Concerning predictive approaches, a large number of methods have been developed or applied to calculate the related Gibbs energy of ion hydration which is easier to quantity than entropy because it connects directly with the equilibrium constant. These include the Born equation, [7][8][9] Kirkwood-Buff theory, 10,[START_REF] Ben-Naim | Statistical Thermodynamics for Chemists and Biochemists[END_REF] perturbation methods, [START_REF] Straatsma | [END_REF][13][14][15][16][17][18][19][20][21] quasichemical theory, [22][23][24] scaled particle theory, 25 and correlation functions. 26,27 Hydration entropy calculations are less commonly reported for ion hydration and may be evaluated from the the Born equation, 8 the temperature derivative of the Gibbs energy change, 21 subtracting off the enthalpy change from the Gibbs energy change, [28][29][30] scaled-particle theory, 25,31 correlation functions, 26,27 or Langevin dipoles. 32,33 We have developed a method that is capable of calculating the entropy of a dilute aqueous solution from computer simulation which has so far been applied to noble gases 34 and small organic solutes. 35 It partitions the entropy into translational and rotational degrees of freedom for every molecule in the system. Configurational space in translation is discretised by the size of the water molecule into distinct minima whose entropy is derived from the force magnitude measured in a computer simulation. Rotational space for water molecules is discretised by hydrogen bonds between the solute and water molecules. Here we extend this approach to calculate the entropy of ion hydration as well as the Gibbs energy and enthalpy. Similar to most computational studies, we examine the alkali metal cations and halide anions because of their simplicity and abundance of computational and experimental data. We examine how the entropy components vary for each ion and use this knowledge to understand commonly used concepts such as kosmotropes and chaotropes and give further insight into the Hofmeister series. 4,[START_REF] Kunz | Specific Ion Effects[END_REF]36 Method Gibbs Energy, Entropy and Enthalpy of Hydration. Using the same theory derived elsewhere for the hydration of noble-gas solutes in water, 34 the standard Gibbs energy, entropy and enthalpy of hydration at temperature T for the transfer of an ion from the gas phase to a dilute aqueous solution at the same concentration are

∆G • X+w = E sim X,w -E sim w(l) -3k B T -k B T ln    1 V w 2k B T F X(aq)
3 Ω or w(X)

Ω or w(l) j=x,y,z

F j w(l) F j w(X) τ j w(l) τ j w(X) N w(X)    (1)
∆H • X+w = E sim X,w -E sim w(l) -k B T (3/2 + α w T ) (2)
∆S • X+w = k B ln    1 V w 2k B T F X(aq) 3 Ω or w(X)
Ω or w(l) j=x,y,z

F j w(l) F j w(X) τ j w(l) τ j w(X) N w(X)    + k B (3/2 -α w T ) (3)
Eq. 2 is exact. Eqs. 1 and 3 assume that the environment of every molecule can be approximated by mean-field harmonic potentials and only include the contribution of water in the first hydration shell of the ions. E sim X,w and E sim w(l) are the energies of the solute in water and of bulk water, k B is Boltzmann's constant, V w = 30.0 Å3 is the volume per water molecule at temperature 298 K and pressure 1 bar, F X(aq) is half the force magnitude on the solute, Ω or w(X) and Ω or w(l)

are the number of orientational minima of a water molecule in the ion's first shell and in bulk, F j w(l) and τ j w(l) are half the force and torque magnitudes on water for the j = x, y, z axes, N w(X)

is the number of water molecules perturbed by the solute, and α w is water's thermal expansion coefficient 37 which is calculated to be 2.8 × 10 -4 K -1 for TIP4P/2005 at 298 K, close to the experimental value 38 of 2.57 × 10 -4 K -1 . We note the presence of the k B T 2 α w term in Eq. 2, amounting to 0.2 kJ mol -1 at 298 K, is usually ignored when enthalpy changes for solvation at constant concentration are evaluated from simulation energies.

As in earlier work, 34,39 Ω or w(l) for bulk water is evaluated using the equation

Ω or w(l) = N w w (N w w -1) 2
N w w -2 N w w 2 (4)
where N w w is the number of first-shell water molecules, defined as being within the first minimum of the oxygen-oxygen radial distribution function g(r OO ) for bulk water which lies at r = 3.4 Å.

This equation can be interpreted as the total number of uncorrelated ways that a water can donate its two hydrogens, dividing by the symmetry number 2, and multiplying by the probability that the neighbouring water has the correct orientation to accept a hydrogen, squared because the donating water has two hydrogen atoms. It represents the number of hydrogenbond arrangements per water molecule. In this work for ions, an extension of the theory to calculate Ω or w(X) is necessary because ions disrupt the neighbouring hydrogen-bond structure much more than noble gases. We again assume that only first-shell water is perturbed by the ion. The number of first-shell water molecules N w(X) is defined as those within the first minimum of the ion-oxygen radial distribution function g(r XO ). The first-shell assumption is based on observations from experiment, 36,[40][41][42] other theoretical treatment, 43 and our own analysis of second-shell water molecules whose entropies we found to be very similar to those in bulk. We did observe more significant perturbations extending into the second-shell for divalent cations but they are not the study of this work. For noble gases, 34 we had observed and consequently 4 assumed that first-shell water molecules have equal probabilities of forming hydrogen bonds to any other water molecule, whether in the first shell or in bulk, and that no hydrogen bonds are formed with the solute. There we defined a hydrogen bond as between a hydrogen and its nearest non-bonded oxygen. For ions, the hydrogen-bond pattern of water would be expected to be very different: first-shell water molecules would have different probabilities to form hydrogen bonds to the ion, to other first-shell water molecules or to bulk water. We denote these probabilities as p i HB where i is the species involved in the hydrogen bond. These probabilities are defined as

p i HB = N i HB N i w(X) (5)
where N i w(X) is the number of water molecules in the first-shell of species i, defined as before, and N i HB is the number that form a hydrogen bond. N i HB may be further decomposed into the number of hydrogen bonds that the water accepts N i a and donates N i d from species i such that

N i HB = N i a + N i d .
However, this decomposition is only carried out for interest and not needed in the entropy calculation. Here we use a more general force-based hydrogen-bond definition than those commonly used based on distance 44,45 or energy 46 because we need to take into account the differing abilities of water and ions to accept, and because these other definitions contain arbitrary cut-off parameters that would give arbitrary entropies. We had employed this force definition in earlier work on the entropies of small organic solutes in water. 35 We define a hydrogen bond between a hydrogen and the atom with which it has the most attractive force.

Given that the hydrogen atoms of water do not have Lennard-Jones radii for the SPC/E force field used here, only the Coulombic term contributes to the force. This amounts to saying that a hydrogen forms a hydrogen bond to the acceptor atom A with the largest value of q A /r 2 AH , where r AH is their distance of separation, and we neglect donation to cations because they have repulsive interactions with the hydrogen atoms of water. We define an effective coordination number N eff w for a first-shell water whereby each value of N i w(X) is weighted by its probability p i HB relative to the largest probability p i,max HB summed over the three types of species i, namely X, first-shell water w(X), and bulk water w(l), giving

N eff w = i=X,w(X),w(l) p i HB p i,max HB × N i w(X) = i=X,w(X),w(l) N i HB p i,max HB (6)
It can be seen that N eff w = N w w when the only species i is bulk water. N eff w is used in Eq. 4 in place of N w w to calculate the number of uncorrelated orientations. We retain the bulk probability that a neighbouring water molecule can accept a hydrogen bond because the effect of the ion on these probabilities has already been taken into account in defining N eff w . The only other change we make to Eq. 4 is in the probability that an ion can accept because the spherical anions considered here, unlike water, can always accept a hydrogen bond. Hence the probability that a neighbouring species can accept is not to the power 2 but rather to the power 2 -p X HB . This gives the final equation

Ω or w(X) = N eff w (N eff w -1) 2
N eff w -2 N eff w 2-p X HB (7)
6 For later analysis, we partition the entropy change in Eq. 3 into four terms ion :

∆S • X = k B ln 1 V w 2k B T F X(aq) 3 + k B (3/2 -α w T ) (8)
water vibrational :

∆S vib w = N w(X) k B ln j=x,y,z F j w(l) F j w(X) (9)
water librational :

∆S lib w = N w(X) k B ln j=x,y,z τ j w(l) τ j w(X) (10)
water orientational :

∆S or w = N w(X) k B ln Ω or w(X) Ω or w(l) . (11)
Molecular Dynamics Simulations. The ions examined in this study are the monovalent cations Li + , Na + , K + , Rb + , and Cs + , and anions F -, Cl -, Br -, and I -. The Lennard-Jones parameters for the ions, listed in Table 1, are the parameters used in the parm99 force field as implemented in AMBER 9. 14,[47][48][49] Simulations of bulk water in a cubic periodic box with 375

SPC/E water 50 and dilute solutions containing in addition one type of each ion are carried out using the AMBER9 software package. 51 500 steps of steepest-descent minimisation is followed by a molecular dynamics simulation in the N V T ensemble at 298 K using a Langevin thermostat. The systems are equilibrated for 2 ns in the N P T ensemble using a Berendsen barostat.

The force and torque magnitudes on the ion and water, coordination numbers N w(X) and N i w(X) , and hydrogen-bond counts N i a and N i d are averaged over a further 1 ns. Longer simulations of 20 ns were used to calculate the energies E sim X,w and E sim w(l) . In evaluating the energies, the size-dependent self-potential 28,[START_REF] Hummer | [END_REF]53 does not need to be accounted for separately 54 because it is already in the particle mesh Ewald implemented in AMBER 9. The thermodynamic finite-size correction 30,[START_REF] Hummer | [END_REF]55

Results

Entropy Change of the Ion. Table 2 contains the ions' force magnitudes F X(aq) and their contribution ∆S • X to the entropies of hydration which are evaluated using Eq. 8. As expected, smaller ions with higher charge densities experience higher force magnitudes and lose more entropy. Anions lose more entropy than their isoelectronic cations because they can get closer to water's hydrogen than cations can to water's oxygen. For similar reasons, an anion's entropy loss is more sensitive to ionic radius (Table 1) than that of cations. To provide more insight into these entropies, we list in Table 2 the free volumes v tr X(aq) = 2k B T /F 3 X(aq) and the ion's contribution to the partial entropy in solution, 34 calculated using s vib X(aq v tr X(aq) values for ions are much smaller than those of aqueous noble gas solutes 34 which range from 1.8 Å3 xenon to 5 Å3 for neon tobecause their force magnitudes can be as low as a third that of water's. The partial entropy of aqueous ions is seen to increase with size because of both the ion's larger v tr X(aq) and its larger mass but this dependence on mass cancels with the gas-phase entropy for solvation.

) = k B ln(v tr X(aq) /Λ 3 X )+3k B +α w k B T where Λ X = h/(2πm X k B T )
Entropy Change of Water. The radial distribution functions g(r OX ) used to calculate N w(X)

are shown in Figs. 1(a) and (b) and the cut-offs r cut XO are listed in Table 3. The resulting values of N w(X) are also listed in Table 3 which are all consistent with experimental values. [START_REF] Ohtaki | [END_REF]18 The next set of figures in Table 3 are the x, y and z components of ∆S vib w and ∆S lib w calculated using Eqs. 9 and 10. Many coordinate frames are possible for water. To be consistent with our approach used for solutes, 35,39 in contrast with the initial formulation, 57, 58 and also adopted by Zielkiewicz, 59 the x, y and z axes are the principal axes ordered from lowest to highest moments of inertia. Thus the x axis is parallel to the hydrogen-hydrogen vector, the y axis is along the dipole moment, and the z axis is perpendicular to water's plane. Unlike the case of the noble-gas solutes for which ∆S vib w and ∆S lib w increased marginally and essentially isotropically, we observe a much richer behaviour for first-shell water around ions. First, anions' first-shell water loses much more vibrational entropy than cations'. Second, anions' first-shell water loses a substantial amount of librational entropy but cations' hardly lose any. Third, the size of the entropy loss is greater for a higher charge density. The first two trends reflect the well-known asymmetry of the water molecule and how cations and anions interact differentially with it. 15,28,[60][61][62][63] An anion interacts with water's hydrogens, pinning the water down translationally and rotationally. A cation interacts more weakly with water's less accessible and more central oxygen, such that the translational and rotational restriction is less. If we consider the directional anisotropy in more detail, the greatest translational loss for water around Li + is along the dipole axis, consistent with a structure having water's dipole pointing away from Li + . This trend weakens for larger ions and for the larger cations we observe a slight increase in all translational components, similar to what was seen for noble gases. 34 This increase is largest in the z component of water,

suggesting that water molecules are partially lying down on the surface of the cations, again, as they would for a non-polar solute. Little anisotropy is seen in the translational entropy loss for anions which indicates that the interaction with the hydrogen can pin the whole molecule down equally in every translational direction. Concerning libration for anions, the restriction is greatest about the y and z axes but less about the x axis because rotation about the hydrogenhydrogen vector has less effect on the position of the hydrogen.

Next we consider the behaviour of the quantities needed to evaluate the orientational entropy. In Table 4 are listed coordination numbers

N X w(X) , N w(X) w(X) , and N w(l)
w(X) of the three types of species X, w(X) and w(l) surrounding a first shell-water. Also included in Table 4 are the number of times these waters accept N i a and donate N i d , as well as the hydrogen-bond probabilities p i HB calculated using Eq. 5. The ion is always in the coordination shell of first-shell water by definition. We observe that N w(X) w(X) decreases with ion size, as does N w(l) w(X) to a smaller extent. This is because of the flatter surface of larger ions that exclude more nearby water molecules from a first-shell water. Concerning hydrogen-bond probabilities, water may donate to anions and other water molecules but not to cations. We note that N F - d is actually calculated to be 1.03, indicating that on a small number of occasions, both of a water's hydrogens point to the F -. We neglect these instances and set N F - d = 1. We next note for first-shell waters that the number of donors equals the number of acceptors by definition. In addition, there is a strikingly small number of hydrogen bonds between first-shell water molecules, a number that increases with solute size. This is because triangles of hydrogen bonds and ion-water interactions are much more strained compared to the more common pentagons, hexagons and heptagons. Non- The ions here lie somewhere between these two extremes depending on their size. To the bulk, first-shell waters donate much more around cations than they do around anions, and conversely for accepting hydrogen bonds.

This effect is consistent with the polarisation of water by ions and again decreases with ion size. Yet even for the largest ions, the ratio between donors and acceptors is still 2:1. Thus the orientational ordering extends beyond the first-shell but we condense this into the first-shell.

Interestingly, water around all anions almost always has one HB to the bulk, this being the hydrogen not interacting with the ion or with first-shell water. In contrast, for cations, the number of first-shell waters donating to bulk waters is nearer two, particularly for the smaller ions. The hydrogen-bond probabilities represent the combination of these effects for donors and acceptors normalised by the coordination numbers. They are close to the probability for bulk water which is p HB = 0.85, calculated using Eq. 5 and the values of N w w = 4.73 from Table 5 and N w HB = 4.

To get the orientational entropy for water around ions, Table 5 contains the effective coordination numbers N eff w calculated using Eq. 6 and the numbers of orientations Ω or w(X) calculated using Eq. 7. Ω or w(X) can be converted into the partial orientational entropy

s or w(X) = k B ln Ω or w(X)
per water molecule which is also in Table 5. The last two columns in Table 5 Entropy, Enthalpy and Gibbs Energy of Hydration. A summary of the four entropy components calculated using Eqs. 8-11 are given in Table 6. These data sum up a key result that the loss of entropy of water around anions is dominated by vibrational and librational entropy while that of water around cations is dominated by orientational entropy. Also in Table 6 is the total ∆S • X+w from Eq. 3 together with ∆S • X+w from exponential averaging for the cations with the same Lennard-Jones parameters in SPC water, 21 and from the Langevin-dipoles method. 32 Exponential averaging measures the free energy difference between the hydrated and gas-phase ions. The Langevin-dipoles method is similar in spirit to ours in that they both evaluate the entropy directly from a single simulation. It is made up of a term dependent on the surface area and a term quantifying rotational restriction of the solvent dipoles. Incidentally, both sets of authors 21,32 compare their entropy changes at constant concentration with those from experiment 64,[START_REF] Burgess | Metal Ions in Solution[END_REF] for the solvation process at 1 bar in the gas-phase to 1 M in aqueous solution for which the entropy change is 26.5 J K the worst is K + which is too positive by 10 J K -1 mol -1 . Similarly, there is good agreement with the Langevin dipoles method for the cations, with all values within 9 J K -1 mol -1 . Concerning the anions, the agreement is worse with the Langevin dipoles. The overall trend is correct but our values are spread over a wider range than those using the Langevin dipoles method, particularly for the larger anions. As will be clearer when comparing with experiment in the next section, it appears that our current model is too generous in the orientational entropy of the anions which are able to accept hydrogen bonds. We cannot quantitatively compare with the study of Horinek et al. 30 because they do not report any values representative of comparable monovalent ions. Moreover, their inclusion of the surface-potential entropy term of 96.5 J K -1 mol -1 appears to lead to unrealistically large entropy losses.

Finally in Table 7 we examine how ∆S probably reflects the experimental data 64,[START_REF] Burgess | Metal Ions in Solution[END_REF] that Aqvist used 14 to derive the cations' Lennard-Jones parameters whose free energy values are about 40 kJ mol -1 more negative than the ones being compared to here. If we also consider the Born correction of 11 kJ mol -1 to account for ion-water interactions beyond the spherical boundary, 14,21 then these terms together account for the discrepancy. For anions, the force field is probably partly responsible, particularly in light of the too-negative enthalpy changes. The disagreement in entropy with values from the Langevin dipoles method suggests that the model for the orientational entropy is also responsible for the difference, probably because an overestimation of the gain in Ω or X(aq) when a hydrogen donates to the ion. Clearly, an improved model for orientational entropy is required, an issue we discuss later.

Discussion

The increasing entropy loss with increasing charge density, either because of a larger charge or a smaller radius, is well known and is the basis of the Born equation. 7 Our results are in accord with this trend. However, an important aspect of ionic hydration that the Born equation cannot explain is the nature of the entropy loss and how it might differ between cations, anions and non-polar solutes. One striking difference between non-polar solutes and ions concerns the heat capacity of solvation. It is generally negative for ions and positive for non-polar solutes. 68,69 This is understood in terms of weaker, more distorted water-water hydrogen bonds in ionic solutions and stronger, more linear water-water hydrogen bonds in non-polar solutions. This is the reverse of the solute-water interactions which are stronger for ions than non-polar solutes. The influence on heat capacity of the water-water interactions, being more numerous than solute-water interactions, is evidently dominant. This different behaviour between ions and non-polar solutes would be a puzzle Lynden-Bell and Rasaiah. 28 The authors calculate the entropy of hydration as a function of ion size and charge by letting the solute charge and radius be dynamic variables. They observed a greater entropy loss for ions of either larger charge or smaller charge, and only a moderate loss for ions of intermediate charges, whether a cation or anion. They also observe maxima in the entropy loss as a function of solute radius for ions but not for non-polar solutes. Similar trends have also been reported by others. 25,26,70 These curves represent a balance between two kinds of entropy loss, non-polar and ionic. Concerning the hydration of cations and anions, there is a well-known difference 15,28,60,62,63,[START_REF] Krestov | Thermodynamics of Solvation: solution and dissolution, ions and solvents, structure and energetics[END_REF][START_REF] Collins | [END_REF] arising from the asymmetric structure of the water molecule. The two positively charged hydrogens are more accessible than the negatively charged oxygen, giving it a positive electrostatic potential at its surface. 15,62 Consequently, anions are better hydrated than cations, meaning that the anion has a more favourable Gibbs energy of solvation than a cation of the same size. Furthermore, there is a weaker dependence of the entropy loss on cation size than anion size. 28,61 Clearly, anions interact more strongly with water's hydrogens than cations do with water's oxygen.

Our results here (Table 6) and elsewhere 34 let us make the following interpretation of the entropy losses: they are a consequence of stronger solute-water hydrogen bonds in the anionic case and of fewer hydrogen-bond arrangements in the non-polar 34 molecules have difficulty hydrogen-bonding with each other and accepting from bulk water, owing to the cation's strong attraction for the oxygens of first-shell waters. This effect for anions is offset by anions' ability to always accept hydrogens which can actually increase the number of hydrogen-bond arrangements, although perhaps excessively in the current model.

An explanation of the entropy loss for cations along similar lines was given by Singh and

Warshel 33 who used the Langevin dipoles method. They found that the entropy loss for Na + , K + and Ca 2+ is in solvent configurational space rather than in solvent vibration modes because water molecules around ions have a narrower variation in their dipole orientations but similar coordinate fluctuations to those in bulk. They did not consider anions. The insensitivity of ∆S or w for cations is also noteworthy. It occurs because the number of first-shell water molecules increases more for cations than for anions. This in turn is because of the greater variation in size in the cations than in the anions (Table 1). Thus, even though the entropy losses per water molecule are reduced for larger ions of either type, the total entropy loss is not as diminished for cations because more water molecules are affected (Table 5). This effect is even stronger for non-polar solutes. 34 They reduce the entropy per water molecule even less but may have up to 21 first-shell water molecules in the case of xenon. This helps explain the non-intuitively greater entropy loss of two argon atoms compared to that of one K + and one Cl -ion. [START_REF] Friedman | Aqueous solutions of simple electrolytes[END_REF][START_REF] Ben-Amotz | [END_REF] In the absence of a full understanding of the entropy of solutions, key concepts used by many to understand ionic hydration are structure-making kosmotropes and structure-breaking chaotropes. [START_REF] Kunz | Specific Ion Effects[END_REF][START_REF] Zhang | [END_REF][START_REF] Collins | [END_REF]75 Care must be taken when connecting them with entropy because the negative entropy of hydration for all ions would suggest that they are all kosmotropes. Rather, one must look at the water's entropy component to make this connection, and there are different ways of partitioning the entropy of solutions between water and the solute. 34,39 Our approach treats all molecules equivalently and assigns entropy to each molecule in the mean-field environment of their neighbours. The solute and solvent entropy are still both negative for cations 16 and anions. However, if one focusses on water's vibrational and librational entropy (Table 6), smaller ions, having negative entropy terms, match the term kosmotropes, and larger ions and non-polar solutes, 34 having positive entropy terms, match the term chaotropes. The effect is stronger for anions than for cations. The trend is consistent with the original naming of structure-makers and structure-breakers by Frank and Evans 56 and the zero-points are consistent with Na + and Cl -being intermediate. 76 An alternative approach is to assign the gas-phase entropy to the solute and an excluded-volume entropy penalty to the solvent to compensate for the solute's greater freedom. 25,27,70,[START_REF] Ben-Amotz | [END_REF] This gives quite a different interpretation of the entropy change of solvation: it arises from an almost complete cancellation between large solute-water and water-water hydration terms. For ions, the solute-water term is negative, indicating that water molecules are very restricted relative to the solute. This is the term that corresponds with kosmotropic and chaotropic behaviour because ions of higher charge density order the surrounding water more relative to the ion. The water-water term, though, is positive because the ion disrupts the water structure relative to how water is normally arranged relative to other water molecules. For non-polar solutes, both these effects are strongly reduced and even have the opposite sign such that the solute-water term is positive and the water-water term negative.

If one assesses water-structuring based on this term, 27, 77-80 contradictory results ensue: ions of higher charge density, namely kosmotropes, break down the structure of water relative to water rather than making it; chaotropes and non-polar solutes increase water-water structure.

Applying the concept of water-structuring beyond how it was originally intended has led many to question the validity of the concepts of kosmotropes and chaotropes. 75,79,80 Nonetheless, the structure-making concept does not tell the full story in relation to the entropy of solutions because it ignores, depending on which entropy decomposition one does, the number of hydrogen-bond arrangements or the relative sizes of the water-water and water-solute terms.

Only a full analysis of the entropy can clarify these issues. 17 Ions are known to have strong effects on solution properties. For example, kosmotropes increase the solution's viscosity and surface tension. [START_REF] Kunz | Specific Ion Effects[END_REF][START_REF] Collins | [END_REF] A particularly interesting and unresolved effect of ions relates to the Hofmeister series [4][START_REF] Kunz | Specific Ion Effects[END_REF][START_REF] Zhang | [END_REF]36 which connects protein stability and solubility in aqueous solution with an approximate sequence of anions and cations. Anionic kosmotropes increase protein stability and promote aggregation, consistent with anions' strong attraction for water. Anionic chaotropes interact more strongly with non-polar and chaotropic parts of the protein such as ammonium and guanidinium groups and thereby solvate it. [START_REF] Collins | [END_REF] Cations, however, have a mixed effect: strong chaotropes enhance protein unfolding and solubility, similar to anionic chaotropes, but so do strong kosmotropes. [START_REF] Kunz | Specific Ion Effects[END_REF] A theory for this in terms of enthalpy is based on the like-attracts-like rule in that a protein's anionic groups are kosmotropes but its cationic groups are chaotropes; 5, 76 strong kosmotropic and strong chaotropic cations interact with the kosmotropic and non-polar parts of the protein, respectively, while intermediate cations interact less strongly. We suspect this difference between anions and cations also has entropic origins in terms of the number of hydrogen-bond arrangements. Anions increase this number while cations, similar to non-polar molecules, decrease it, meaning that cations may favour the solvation of non-polar surfaces of the protein. In other words, the surface of the first hydration shell of kosmotropic cations appears mostly hydrogen-like, giving it a more hydrophobic guise, similar to strong chaotropes such as guanidinium and urea. Thus kosmotropic and chaotropic cations would both have reasons to unfold or solvate proteins, thereby explaining the mixed order of cations with respect to kosmotropic and chaotropic behaviour in the Hofmeister series and the weaker influence of cations than anions. It is difficult to resolve unambiguously the myriad of competing enthalpic and entropic effects involving ions, water and protein atoms without studying proteins in aqueous electrolytes directly.

More accurate and extensive studies of these effects at the level of individual ions can be made using better force fields and models for the entropy of aqueous solutions. A number of force fields for ions have recently been developed, of both the fixed-charged 18, 66 and polarisable nature. 81 As we noted earlier, the inclusion of quantum effects would be necessary for ions of higher charge density, particularly for transition-metal complexes with more covalent interactions. Regarding methods to calculate the entropy of aqueous solutions, one of the deficiencies in our theory is the dependence on arbitrary cut-offs in defining hydration shells. This is not as serious a problem as fitted parameters which inevitably would give better agreement, but are still unsatisfactory. These also make more problematic the study of aspherical solutes and longer-range entropic contributions. Another deficiency is the restriction to one type of water environment which had been based on the prevailing view that water structure is essentially tetrahedral. 4,82 This leads to a more awkward mean-field theory for the orientational entropy that clearly averages too much over variability in water structure around ions. Recent work by us has found a way past both these difficulties. We have proposed a new hydrogen-bond definition that has no arbitrary parameters: a hydrogen donates to its nearest oxygen as long as that oxygen-hydrogen distance is the shortest between the two water molecules involved. 83 This definition accurately resolves the hydrogen-bonding pattern in water and reveals multiple water environments according to the number of hydrogen bonds formed rather than a single distorted tetrahedral arrangement. It can be generalized in the case of different donors and acceptors by using the electrostatic force as we have already done here. It points the way to a much more flexible theory for solutions as a mixture of multiple water and solute environments based on their hydrogen-bond patterns. The development of theory to do this is currently in progress.

Conclusion

A method has been presented to calculate the entropy of hydration for alkali metal cations and halide anions from molecular dynamics simulations of each ion in water and of bulk water. En- tropy changes are derived from the force magnitude on the ion, the force and torque magnitudes and coordination and hydrogen-bond numbers of its first hydration and bulk water molecules.

The entropy results are in reasonable agreement with perturbation calculations for cations and partial agreement with the Langevin dipoles calculations for both types of ions. Agreement in the entropy change with experiment is better for smaller ions but is too negative for larger cations and not negative enough for larger anions. There is a much greater vibrational and librational entropy loss for the anions than cations which is consistent with stronger water-ion interactions and more favourable hydration for anions. The number of orientational minima, or equivalently, the number of hydrogen bond arrangements decreases substantially for cations but increases moderately for anions. This is because cations make it difficult for first-shell water molecules to accept hydrogen-bonds from other water molecules whereas anions affect this little and moreover can always accept a hydrogen. The terms of kosmotropes and chaotropes are shown to apply to the first-shell water's vibrational and librational entropy but not to the overall solution entropy. Finally, the number of hydrogen bond arrangements is proposed to make an important contribution to the Hofmeister ordering of ions in that anions increase this number but cations decrease it, and may help explain kosmotropic cations solvate and unfold proteins and why cations have a weaker influence than anions. "In this paper the hydration properties of several monovalent cations and anions is evaluated using an approximate theoretical expression. I think this paper can be published in Molecular Physics, although I must confess that I have certain doubts about the real practical use of the approach used by the authors."

F

In relation to the referee's doubts about the real practical use of our approach, we are not clear why the referee makes this judgment. The manuscript already demonstrates the practical use of the approach. Granted it is approximate (harmonic approximation, and only two water environments considered, namely first shell, each with a mean-field number of orientations)

and it has only been implemented for a certain class of system (i.e. radially symmetric, dilute ions). However, a person wishing to study solvation has few alternative methods available to understand the entropy of solutions and even these have their own approximations. For example, the method of distribution functions requires approximations such as the omission of water-water correlations to keep dimensionalities at a practical. Despite these limitations, our method gives new and quantitative insights into ionic solvation. Furthermore, it could be directly extended to study the hydration of a solute of arbitrary complexity, again with these same limitations. At the end of the manuscript we allude to a more general approach we are developing using our new hydrogen bond definition that examines multiple water environments and avoids the use of arbitrary cut-off parameters. I also made a brief mention in my seminar at the Liblice Conference that we are not limited to the harmonic approximation. Work addressing these two points is still being carried out and is not yet ready for publication.

"1. The authors should state that Eq.( 1) is not an exact expression but rather an approximation.

In fact the approximate expression divides the free energy into a contribution of the ions and a contribution of the water."

We have made it clear in the manmuscript that Eq. 1 is an approximation. After these equations, we write "Eq. 3 is exact. Eqs. 1 and 2 assume that the environment of every molecule can be approximated by mean-field harmonic potentials and only include the contribution of water in the first hydration shell of the ions."

"2. The authors use the experimental value of α in Eq.( 3) and not the value of α of the water model used by them."

Yes, to be consistent we now use the simulation value of α in Eq. ( 3) rather than the experimental value (0.000257 K -1 ). This has been calculated by the developers of TIP4P/2005 to be 0.000028 K -1 (Abascal and Vega, JCP, 2005, 123, 234505) which is close to the experimental value. This makes a negligible difference to the energy (0.02 kJ mol -1 ) and entropy (0.06 J K -1 mol -1 ). Some entropy values in Table 2 (and the matching second column of Table 6) decreased by 0.1 J/K/mol. No change occurred in the final entropy or enthalpy values, reported to lower precision. We also corrected an error in the size of the αk B T 2 term given in the text to now read 0.2 kJ mol -1 instead of 0.6 kJ mol -1 which was the value of αk B T but with the units of energy. This was a typo and has no effect on the any of the other calculated values.

"3. I do not understand why the mass is included in Table 2. Although the mass of a molecule appears in the expression of the free energy, it is absolutely irrelevant when it comes to describing transfer properties (basically the mass contribution is the same in both phases). This point should be clarified in the text." Yes, to calculate the solvation terms, we do not need the mass. The mass is only needed to calculate the ion-entropy components in solution which is listed in the second-last column of Table 2. However, we agree that the masses of monatomic ions are commonly available data and do not need to be included here. Thus we have removed the masses from Table 2.

"4. The authors use a force field to perform simulations of the ions in water. After using this force field the authors extract the required information to implement their theory (mean forces, torques ...) . Is the force field used by the authors able to describe the hydration free energies and entropies of the ions described in the paper? Could the authors comment on this?"

We compare our entropy results with those calculated elsewhere using thermodynamic integration with the same ion force field for cations in SPC water (Carlsson and Aqvist, ref 15) and obtain good agreement. We have not found equivalent data for anions and nor have we done the calculations using thermodynamic integration. However, we can assess the performance of the ion force field compared to experiment by looking at the enthalpy: the solvation enthalpy values are not negative enough for cations and too negative for anions. These discrepancies would not affect entropies as much because entropy is much less sensitive to force field than entropy e.g. going from Li + to Na + entails a ∆H of 117 kJ mol -1 but a T ∆S of only 7 kJ mol -1 .

Differences in force field parameters for a given ion are much subtler. Entropy is less sensitive because it depends on the logarithm of configurational volume and the associated force field parameters. Enthalpy, on the other hand, has an approximate linear dependence on force field energy and charge parameters.

"5. The agreement between theory and experiment shown in Table 7 is rather impressive.

The authors used theory. I have the feeling that using SPC/E and the force field for the ions described here the hydration properties would not match the experimental results. My feeling is that if you use theory to describe experimental properties using information from a force field that does not reproduce experimental properties, then you are basically fitting the theory parameters to reproduce the experimental results. Could the authors comment on that?"

We do not fit any parameters to reproduce experimental results. We emphasise that all quantities in our entropy calculations are derived from the computer simulation which in turn is only based on force-field parameters. This approach is consistent with our objective of predicting free energy from a simulation without recourse to any fitted parameters. In our theory, there are two arbitrarily assigned parameters which are the cut-off radii at the first minima in the oxygen-oxgyen and oxygen-ion radial distribution functions. These were chosen a priori and not fitted.

Finally, we have rewritten the discussion to make it clearer and better organised and improved the wording throughout the manuscript.

We hope these alterations satisfy you and the referee.

  and cationic cases. In other words, the potential energy surface of aqueous anions have deeper vibrational and librational energy wells while aqueous cations and non-polar solutes have fewer energy wells. This explanation is consistent with the more favourable hydration of anions and the stronger dependence on charge density. The small change in vibrational and librational entropy of cations is consistent with their weaker interaction with water. The number of hydrogen-bond arrangements is reduced for non-polar solutes and cations because of their inability to accept hydrogen bonds from water. The number is lower for cations than non-polar solutes because first-shell water 15

Page 17 of 35 URL

 35

Figure 1 .

 1 Figure 1. Ion-oxygen radial distribution functions g(r XO ) of (a) the five cations and (b) the four anions. Each vertical dashed line is the cut-off r cut XO placed at the first minimum.

Page 6 of 35 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

  is only a few kJ mol -1 and we omit it here. Simulations have SHAKE on all

	bonds and a 2 fs timestep, periodic boundary conditions, 8 Å cutoff and particle mesh Ewald
	with default AMBER parameters. Equilibrated box sizes have sides of ∼22.5 Å.
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	7

  polar solutes, on the other hand, are encaged by and separate from water's hydrogen-bond network, making a intra-first-shell hydrogen bond quite normal.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	10

  N eff w values around ions lie below that of water, primarily because it is difficult for first-shell water molecules to form hydrogen bonds with other first-shell water molecules while simultaneously interacting with the ion. This effect is strongest for the cations and smaller ions. ways dependent on neighbouring water molecules having the correct orientation to be accepted.Consequently, there is a loss of orientational entropy per water molecule around cations but a gain per molecule around anions. The total losses of orientational entropy for cations are similar because the smaller loss per water molecule is offset by the larger coordination number N w(X) . The losses become less negative for larger cations and more positive for larger anions. This is because first-shell waters around larger ions form more hydrogen bonds with each other, an effect in turn made easier by the weaker interaction with the ion and the ion's flatter curvature which places other first-shell water molecules at an angle more compatible with water's

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	When converting these values to Ω or w(X) , the first-shell waters of anions have much larger values
	than those of the cations and even of bulk water because waters that donate to the anions are
	11

are the change in orientational entropy per first-shell water molecule and the total for each ion, each calculated using Eq. 11.

always accepted. Waters around cations or in bulk water do not have this feature and are al-tetrahedral angle of hydrogen-bonding. However, this gain in orientational entropy for anions may be exaggerated as shown by a comparison with other methods to calculate entropy and with experiment which we next examine.

  -1 mol -1 more negative. Agreement with values obtained from exponential averaging 21 is good for the cations. Most values are within 6 J K -1 mol -1 and

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	12

  20 They are listed for the same solvation process at constant concentration in the gas and liquid phases and include enthalpy and entropy changes. The general trends versus experiment are correct in that small ions have negative ∆S • X+w , ∆H • X+w and ∆G • X+w which become less negative for larger ions. However, there are some discrepancies. The agreement with the entropies becomes worse for larger ions, not being negative enough for larger cations and too negative for larger anions. Given the agreement for cations with exponential averaging21 described earlier, the difference for cations is likely to be a force field issue. The enthalpy changes of cations match the trend of experiment but are consistently 25 kJ mol -1 too positive. This

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	13

• X+w together with ∆H • X+w from Eq. 2 and ∆G • X+w from Eq. 1 all compare with experiment as a test of the whole modelling procedure. Following the rationale of Joung and Cheetham, 66 we use the experimental values of Schmid. 67 These values do not include a contribution from the air-water surface potential which is consistent with values calculated for infinitely periodic systems using Ewald with tin-foil boundary conditions such as the one here. Only in the case of a spherical boundary potential are such terms required. 19 In addition, the values of other works are very close to Schmid et al.'s within an additive offset.

  of entropy loss. The nature of the differing entropy losses is evident in the insightful study of

	Molecular Physics	Page 14 of 35
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	3, 69 if there was only one kind
	14

URL: http://mc.manuscriptcentral.com/tandf/tmph

19 Page 19 of 35 URL

	Molecular Physics
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

: http://mc.manuscriptcentral.com/tandf/tmph

Table 1 .

 1 Lennard Jones Parameters of the Ions cation r min /2 ( Å) ǫ (kcal mol -1 ) anion r min /2 ( Å) ǫ (kcal mol -1 ) Values in parentheses are entropies for the quantum harmonic oscillator57,58

	Li +	1.137	0.01830
	Na +	1.868	0.00277	F -	1.75	0.061
	K +	2.658	0.00033	Cl -	1.95	0.265
	Rb +	2.956	0.00017	Br -	2.22	0.320
	Cs +	3.395	0.00008	I -	2.35	0.400
	Table 2. Average Force Magnitude and Entropy Components of the Ions
	ion	F X(aq) (10 -10 N) (10 -30 Å3 ) (J K -1 mol -1 ) (J K -1 mol -1 ) v tr X(aq) s vib, a ∆S • X X(aq)
	Li +	2.80	0.03	17.6 (19.3)	-47.0
	Na +	2.18	0.05	38.8 (39.1)	-40.8
	K +	1.75	0.10	50.9 (51.0)	-35.3
	Rb +	1.61	0.13	62.7 (62.8)	-33.2
	Cs +	1.41	0.20	71.5 (71.5)	-30.0
	F -	3.09	0.02	27.7 (28.5)	-49.5
	Cl -	2.36	0.04	42.2 (42.5)	-42.7
	Br -	1.97	0.07	56.8 (56.9)	-38.3
	I -	1.83	0.09	64.4 (64.5)	-36.5

a

Table 3 .

 3 Cut-off Distance, Ion Coordination Number and Change in Vibrational and Librational

	Entropy of First-Shell Water
	ion	r cut XO ( Å)	N w(X) ∆S vib,x w (J K -1 mol -1 ) S vib,y w S vib,z w	∆S lib,x w (J K -1 mol -1 ) S lib,y w S lib,z w
	Li +	2.9	4.4	-2.7	-11.1 -2.0	-0.9	0.7 -1.7
	Na + 3.3	5.8	1.0	-1.8	1.0	0.0	0.8	0.5
	K +	3.6	6.7	2.0	2.1	3.7	0.1	0.2	0.4
	Rb + 3.8	7.4	1.7	2.7	4.2	0.3	0.1	0.5
	Cs + 4.0	8.3	1.7	3.3	5.1	-0.6	-0.8 -0.4
	F -	3.2	6.1	-10.8	-8.0 -9.0	-10.6 -15.3 -14.0
	Cl -3.6	6.9	-2.4	-2.3 -2.7	-3.9	-5.6 -4.9
	Br -3.9	7.5	1.7	0.4	0.5	-0.1	0.4	0.3
	I -	4.0	7.6	3.2	1.3	1.2	1.3	2.6	2.1
							27

Table 4 .

 4 Coordination and Hydrogen-Bond Numbers and Probabilities for First-Shell Water

	ion	N X w(X)	N	w(X) w(X)	N	w(l) w(X)	N X a	N	w(X) a	N	w(l) a	N X d	N	w(X) d	N	w(l) d	p X HB	p	w(X) HB	p	w(l) HB
	Li +	1	2.38	3.04	0	0.02	0.36	0	0.02	1.98	0	0.02 0.77
	Na +	1	1.57	3.13	0	0.10	0.61	0	0.10	1.90	0	0.12 0.80
	K +	1	1.11	3.05	0	0.23	0.73	0	0.23	1.77	0	0.41 0.82
	Rb +	1	1.12	2.96	0	0.30	0.76	0	0.30	1.70	0	0.54 0.83
	Cs +	1	1.15	2.89	0	0.40	0.80	0	0.40	1.60	0	0.69 0.83
	F -	1	1.43	3.38	0	0.04	1.93 1.00 0.04	0.94 1.00 0.05 0.85
	Cl -	1	0.85	3.26	0	0.08	1.88 0.97 0.08	0.95 0.97 0.19 0.87
	Br -	1	0.65	3.18	0	0.13	1.83 0.92 0.13	0.95 0.92 0.40 0.88
	I -	1	0.56	3.20	0	0.14	1.84 0.89 0.14	0.97 0.89 0.51 0.88

Table 5 .

 5 Effective Coordination Number, Number of Orientations, and Orientational Entropy

	for First-Shell Water
	ion	N eff w	Ω or X(aq)	s or w(X)	∆S or w /N w(X) ∆S or w
					(J K -1 mol -1 )
	Li + 3.10 1.08	0.7	-8.3	-36.2
	Na + 3.37 1.33	2.3	-6.6	-38.0
	K +	3.60 1.56	3.7	-5.3	-35.4
	Rb + 3.69 1.65	4.2	-4.8	-35.3
	Cs + 3.84 1.82	5.0	-4.0	-33.0
	F -	3.86 3.24	9.8	0.8	5.0
	Cl -4.08 3.56 10.6	1.6	11.0
	Br -4.32 3.95 11.4	2.5	18.4
	I -	4.48 4.23 12.0	3.0	23.1
	bulk 4.73 2.94	9.0	0	0
						28

Table 6 .

 6 Hydration Entropy Components and Comparison with Other Simulations (J K -1

	mol -1 )
	ion	∆S • X	∆S vib w	∆S lib w	∆S or w		∆S • X+w
						this work ref 15 ref 16
	Li + -47.0 -15.8 -1.9 -36.2	-101	-104	-93
	Na + -40.8 0.2	1.3 -38.0	-77	-87	-82
	K +	-35.3 7.8	0.8 -35.4	-62	-56	-71
	Rb + -33.2 8.6	0.9 -35.3	-59	-52	-67
	Cs + -30.0 10.1 -1.8 -33.0	-55	-53	-60
	F -	-49.5 -27.9 -40.0	5.0	-112	-	-85
	Cl --42.7 -7.4 -14.4	11.0	-53	-	-70
	Br --38.3 2.6	0.6	18.4	-17	-	-64
	I -	-36.5 5.7	6.1	23.1	-2	-	-58

Table 7 .

 7 Entropy, Enthalpy and Gibbs Energy of Hydration versus Experiment

	Page 29 of 35				Molecular Physics
	ion	∆S • X+w	∆H • X+w	∆G • X+w
		(J K -1 mol -1 )	(kJ mol -1 )	(kJ mol -1 )
		this work ref 55 this work ref 55 this work ref 55
	Li +	-100	-100	-481	-506	-451	-476
	Na +	-77	-68	-364	-391	-340	-371
	K +	-61	-32	-280	-308	-261	-298
	Rb +	-58	-21	-258	-283	-241	-276
	Cs +	-54	-16	-224	-258	-207	-253
	F -	-112	-127	-639	-539	-605	-501
	Cl -	-53	-75	-499	-392	-483	-373
	Br -	-16	-49	-419	-361	-414	-346
	I -	-1	-31	-385	-321	-384	-311
					29

URL: http://mc.manuscriptcentral.com/tandf/tmph Manuscript Title: Prediction and interpretation of the hydration entropies of monovalent cations and anions (Special Issue Paper for the Eight Liblice Conference on the Statistical Mechanics of Liquids). Manuscript Authors: Sheeba Jem Irudayam and Richard H. Henchman Dear Prof. Jackson, Thank you for examining our manuscript and passing on the referee's comments to us. Our responses to the referee are the following:

Acknowledgments

RHH and SJI are funded by EPSRC grant EP/E026222/1 and SJI is funded by an Overseas Research Scholarship and the School of Chemistry at the University of Manchester. 20


MD Stefan Vordenbäumen
email: stefan.vordenbaeumen@med.uni-duesseldorf.de

Daniel Timm

Ellen Bleck

MD Jutta Richter

MD Rebecca Fischer-Betz

MD Gamal Chehab

MD Oliver Sander

MD Matthias Schneider

Altered serum levels of human neutrophil peptides (HNP) and human beta-defensin 2 (hBD2) in Wegener's Granulomatosis

Keywords: human neutrophil peptide, innate immunity, alpha-defensin, betadefensin, antimicrobial peptides, Wegener Granulomatosis

Defensins are highly conserved peptides with antimicrobial and immunomodulatory functions. Due to their chemotactic properties on mononuclear cells, including dendritic cells and macrophages, defensins may contribute to granuloma formation in Wegener's granulomatosis (WG). In order to explore if these peptides might be involved in WG pathogenesis, sera of patients were screened to detect altered defensin-levels. For this purpose, serum and EDTA-blood of patients with WG (n = 17; aged 54.8 ± 15.5 years) and age-and sex-matched healthy controls (n = 24; aged 55.5 ± 16.8 years) were collected. Levels of neutrophil alpha-defensins (human neutrophi peptides, HNP) and beta-defensin (hBD) 2 and 3 in serum were measured by ELISA. By this means, WG patients displayed higher serum levels of hBD2 and HNP when compared to controls. Furthermore, serum hBD2 was raised in patients with meningeal granulomas (n = 4) or in those undergoing treatment with cyclophosphamide (n = 4). In order to detect if increased gene-expression in polymorphonuclear cells is responsible for the elevated defensin-levels, real-time polymerase-chain-reaction with gene-specific primers was performed. Expression of specific mRNA in polymorphonuclear cells was observed for HNP only, but did not parallel HNP serum levels, suggesting that degranulation rather than increased geneexpression may be responsible for increased HNP-serum levels in WG. In conclusion, elevated serum levels of HNP and hBD2 in WG patients suggest a role for both defensins in WG pathogenesis.

Introduction

Defensins initially attracted interest due to their potent antimicrobial functions and highly conserved structure [START_REF] Schneider | Human defensins[END_REF]. Consequently, several immunomodulatory functions like chemotactic effects and the ability to induce the production of cytokines have been identified [START_REF] Yang | Beta-defensins: Linking innate and adaptive immunity through dendritic and t cell ccr6[END_REF][START_REF] Yang | Multiple roles of antimicrobial defensins, cathelicidins, and eosionphil-derived neurotoxin in host defence[END_REF], suggesting that these peptides connect innate and adaptive immune-responses [START_REF] Yang | Beta-defensins: Linking innate and adaptive immunity through dendritic and t cell ccr6[END_REF]. Their importance for autoimmune disorders is supported by the recent finding of modulated defensin expression in patients with systemic lupus erythematosus [START_REF] Sthoeger | High alphadefensin levels in patients with systemic lupus erythematosus[END_REF][START_REF] Vordenbäumen | Elevated levels of human beta-defensin 2 (hbd2) and human neutrophil peptides (hnp) in systemic lupus erythematosus[END_REF]. Human neutrophil peptides (HNP) stem from the neutrophil granula of leukocytes [START_REF] Rehaume | Neutrophil-derived defensins as modulators of innate immune function[END_REF], while other structurally similar alpha-defensins are mainly expressed in the small intestine [START_REF] Mallow | Human enteric defensins[END_REF][START_REF] Vordenbäumen | Defensins are differentially expressed with respect to the anatomic region in the upper gastrointestinal tract of children[END_REF]. Amongst the beta-defensins (hBD), distinguished from alpha-defensins by diverse disulfid bonding, hBD 2 and 3 are upregulated in response to inflammatory stimuli [START_REF] Becker | Cd14-dependent lipopolysaccharide-induced beta-defensin-2 expression in human tracheobronchial epithelium[END_REF][START_REF] Fahlgren | Betadefensin-3 and -4 in intestinal epithelial cells display increased mrna expression in ulcerative colitis[END_REF][START_REF] Vordenbäumen | Defensin-mrna-expression in the upper gastrointestinal tract is modulated in celiac disease and helicobacter-pylori -positive children[END_REF] and are secreted by some epithelial surfaces and blood cells [START_REF] Zhang | Contribution of human alphadefensins 1, 2, and 3 to the anti-hiv-activity of cd8 antiviral factor[END_REF][START_REF] Duits | Expression of beta-defensin 1 and 2 mrna by human monocytes, macrophages and dendritic cells[END_REF].

Wegener granulomatosis (WG) is characterized by a granulomatous inflammation with various degrees of vasculitis leading to tissue necrosis with a predilection for the respiratory tract and the kidneys [START_REF] Moosig | Wegener's granulomatosis: The current view[END_REF]. Antineutrophil cytoplasmatic antibodies (ANCA) directed against proteinase 3 are strongly associated with WG [START_REF] Hagen | Diagnostic value of standardized assays for anti-neutrophil cytoplasmic antibodies in idiopathic systemic vasculitis[END_REF] and have been implicated in the pathogenesis of the disease [START_REF] Csernok | Pathways to anca production: From differentiation of dendritic cells by proteinase 3 to b lymphocyte maturation in wegener's granuloma[END_REF]. They activate neutrophil granulocytes and induce the release of HNP from neutrophils [START_REF] Tanaka | Fc(gamma)riiib allele-sensitive release of alpha-defensins: Anti-neutrophile cytoplasmic antibody-induced release of chemotaxins[END_REF]. By the ability to chemoattract monocytes, macrophages, dendritic cells, and T cells [START_REF] Rehaume | Neutrophil-derived defensins as modulators of innate immune function[END_REF][START_REF] Röhrl | Human b-defensin 2 and 3 and their mouse orthologs induce chemotaxis through interaction with ccr2[END_REF], defensins may contribute to granuloma formation in WG [START_REF] Tanaka | Fc(gamma)riiib allele-sensitive release of alpha-defensins: Anti-neutrophile cytoplasmic antibody-induced release of chemotaxins[END_REF]. However, reports about defensinserum levels in WG patients are lacking. Therefore, our objective was to evaluate the serum level and expression of HNP, hBD2, and hBD3 in WG patients to further explore if defensins might play a role in WG pathogenesis.

Materials and Methods

Patients and samples

Patients were diagnosed according to ACR criteria [START_REF] Leavitt | The american college of rheumatology 1990 criteria for the classification of wegner's granulomatosis[END_REF]. In all consecutive patients, total blood count, creatinine, c-reactive protein (CRP), and ANCA-staining were routinely performed; current and past medication as well as disease manifestations were recorded. Patients were excluded if the presence of infection could not be excluded based on history, physical examination, and urine analysis. Disease activity was determined by the Birmingham Vasculitis Activity Score (BVAS). For performance of defensin-ELISA, 10 ml of Serum was collected from patients and healthy volunteers and frozen to -20°C prior to the assay. For PCR-analysis, 10 ml EDTA blood was obtained and processed within 2 hours. All patients and volunteers gave their full consent. The study was approved by the local ethics committee.

Enzyme-linked immunosorbent assay (ELISA)

Human HNP1-3 ELISA Test Kit (Hycult biotechnology, Uden, The Netherlands) and human -Defensin 2 and 3 ELISA Kit (Phoenix Pharmaceuticals, Burlingame, California, USA) were applied for measurement of HNP and hBD2 and 3 respectively according to the manufacturer's instructions. Briefly, serum was diluted 1:2000 in dilution buffer for HNP, 1:100 fo hBD2 and undiluted for hBD3 prior to transferring duplicate samples to microtitre plates and initiation of recommended incubation procedures. Absorbance was measured at 450nm using the Anthos 2001 ELISAreader (anthos Mikrosysteme, Krefeld, Germany).

Semi-quantitative real-time Polymerase chain reaction (rtPCR)

Polymorphonuclear cells (PMN) were isolated from EDTA-blood using Polymorphprep (Axis-Shield, Oslo, Norway), mRNA was isolated with Rneasy Mini Kit (Qiagen, Hilden, Germany), and reverse transcription was performed using QantiTect Reverse Transcription (Qiagen) according to the manufacturer's

Statistical analysis

Data comparison was carried out using Mann-Whitney-U test. For comparison of continuous variables, two-tailed bivariate correlations and Pearson's coefficient were calculated. P-values < 0.05 were considered significant.

Results

Patient

Defensin measurements

In order to determine serum levels of defensins, ELISA for hBD2, hBD3 and HNP were performed on samples of patients and healthy controls. HBD2 and HNP, but not hBD3, were elevated in patients with WG as opposed to healthy controls (figure 1).

To determine if defensin-levels varied with respect to the disease state, serum levels were compared to clinical data and laboratory values. Thereby, hBD2 was found to be higher in patients with meningeal involvement (n = 4; figure 2), evidenced by the presence of proliferations on magnetic resonance imaging. Furthermore, patients undergoing treatment with cyclophosphamide displayed superior levels of serum HNP (n = 4; figure 2). Next, we assessed whether PMN could be the source of defensins in WG patients. For this purpose, we performed real-time PCR with gene specific primers for the defensins studied. Only HNP-mRNA was detectable, but did not parallel serum HNP measurements. Values between patients and controls were rather equal.

Discussion

ANCA-mediated release of -defensins might contribute to granuloma-formation due to chemotactic properties on monocytes, macrophages, dendritic cells, and T cells [START_REF] Rehaume | Neutrophil-derived defensins as modulators of innate immune function[END_REF][START_REF] Tanaka | Fc(gamma)riiib allele-sensitive release of alpha-defensins: Anti-neutrophile cytoplasmic antibody-induced release of chemotaxins[END_REF]. Similarly, hBD2 and 3 chemoattract macrophages and monocytes [START_REF] Röhrl | Human b-defensin 2 and 3 and their mouse orthologs induce chemotaxis through interaction with ccr2[END_REF].

Interestingly, high -defensin-levels have been reported in the bronchoalveolar lavage fluid and plasma of one WG patient previously [START_REF] Mukae | Neutrophil-related cytokines and neutrophil products in bronchoalveolar lavage fluid[END_REF]. Furthermore, other vasculitic diseases, like microscopic polyangiitis [START_REF] Tamiya | Defensin-and cathepsin g-anca in systemic lupus erythematosus[END_REF] and cutaneous vasculitis [START_REF] Srivastava | Immunology of cutaneous vasculitis associated with both etanercept and infliximab[END_REF] were associated to increased -defensins and hBD2 expression respectively. The present study demonstrates elevated serum values of HNP and hBD2 in WG patients when compared to healthy controls. Meningeal granulomas and cyclophosphamid treatment, which might be found more frequently in sicker patients, were associated with modulated hBD2 and HNP serum levels respectively. Although hampered by the small patient number, these findings support the notion that defensins might be involved in WG pathogenesis and are an interesting target for future research as novel therapeutic options may arise.

Neutrophile granula of leukocytes are the main source of human neutrophile peptides (-defensin 1-4) [START_REF] Rehaume | Neutrophil-derived defensins as modulators of innate immune function[END_REF]. Interestingly, no correlation of HNP serum levels to the leukocyte count or to HNP-mRNA copy-numbers were found. Thus, we speculate that degranulation of neutrophile granulocytes rather than increased neutrophile numbers or increased HNP-gene-expression is responsible for elevated -defensin levels in WG patients. Beta-defensins are expressed in various tissues which might contribute to beta-defensin serum levels [START_REF] Schneider | Human defensins[END_REF]. Based on our observations, PMN are unlikely to be responsible for beta-defensin production, because no specific mRNA was found.

Although the study served to address our objectives, the small patient number and relatively high standard deviations may be regarded as drawbacks to this study.

Conclusions

To the best of our knowledge, this is the first study to systematically evaluate serum defensin-levels in WG. Our results further implicate a role for HNP and hBD2 in WG pathogenesis. Additional studies in this interesting field are needed.

Legends for illustrations

  instructions. Polymerase-chain-reaction (PCR) with real time measurement of fluorescence was carried out on the 7300 Real-time PCR system (Applied Biosystems, Foster City, CA, USA) with 0.3 µM gene-specific, exon-spanning primers for SYBR green PCR Kit (Qiagen). Results were relatively quantified using glyceraldehyde 3 phosphate dehydrogenase (fw: 3'-CCAGCCGAGCCACATCGCTC-5', rv: 5'-ATGAGCCCCAGCCTTCTCCAT-3') as internal and reference RNA (Stratagene, La Jolla, CA, USA) as external standard according to the -ΔΔCT-method.

Fig. 1

 1 Fig.1Serum levels of human -defensin 2 (hBD2) and human neutrophil peptides (HNP) in patients with Wegener's granulomatosis (WG; n = 17) and healthy controls (n= 24). Significant differences as determined by Mann-Whitney-U test (p<.05) are indicated (*).

Fig. 2 a

 2 Fig. 2 a) Serum levels of human -defensin 2 (hBD2) in patients with (n = 4) or without (n = 13) meningeal proliferations demonstrated by magnetic resonance imaging in patients with Wegener's granulomatosis. b) Serum levels of human neutrophil peptides (HNP) in patients undergoing treatment with cyclophosphamide (CYC; n = 4) compared to those without (no CYC). Significant results as determined by Mann-Whitney-U test are indicated.

Table 1 -Patient characteristics

 1

	Patient characteristics
	number	17
	age (years ± SD)	54.8 ± 15.5
	%female	47.1
	Laboratory values (mean ± SD)
	RBC (/µl)	4,3 ± 0,7
	WBC (/µl)	8.5 ± 4.4
	thrombocytes/µl	278,471 ± 101
	CRP (mg/dl)	1.6 ± 1.9
	creatinine (mg/dl)	1.5 ± 1.0
	PR3-ELISA (U/l)	90.36 ± 100
	Disease activity parameters (mean ± SD)
	BVAS	12.6 ± 14.6
	Organ involvement (n (%))
	renal	8 (47.1)
	upper airway	9 (52.9)
	lower airway	5 (29.4)
	eyes	1 (5.9)
	meningeal	4 (23.5)	p<.05 serum hBD2
	peripheral nerves	2 (11.8)
	Treatment (n (%))
	glucocorticoids	11 (64.7)
	methotrexate	3 (17.6)
	azathioprine	7 (41.2)
	anti-TNF	2 (11.8)
	rituximab	1 (5.9)
	cyclophosphamide	4 (23.5)	p<.05 serum HNP
	Characteristics, laboratory values, disease extend and treatment of patients with
	Wegener granulomatosis. Organ involvement was defined according to ACR criteria.
	RBC: red blood cell count; WBC: white blood cell count; CRP: c-reactive-protein.
	PR3-ELISA: Proteinase3 Enzyme-linked immunosorbent assay; BVAS: Birmingham
	Vasculitis Activitiy Score; meningeal: meningeal involvement as defined by magnetic
	resonance imaging; significant differences determined by Mann-Whitney-U test are
	marked.

financial support information: this work was supported by an unconditional grant form the "Hiller-Stiftung" word count: 1151


Abderrahman Iggidr
email: abderrahman.iggidr@inria.fr

Gauthier Sallet
email: gauthier.sallet@inria.frbergetsanou@yahoo.fr

Berge Tsanou

Global stability analysis of a metapopulation SIS epidemic model *

Keywords: Metapopulation models, SIS models, nonlinear dynamical systems, global stability, monotone systems

 that a SIS type model in a moverstayer epidemic model is globally asymptotically stable is confirmed analytically. If the basic reproduction number R 0 ≤ 1, then the disease free equilibrium is globally asymptotically stable. If R 0 > 1, then there exists a unique endemic equilibrium which is globally asymptotically stable on the nonnegative orthant minus the stable manifold of the disease free equilibrium.

1 Introduction [START_REF] Arino | A multi-city epidemic model[END_REF] considered the so-called "mover-stayer" model adapted from [START_REF] Sattenspiel | A structured epidemic model incorporating geographic mobility among regions[END_REF] by adding the age structure of a stationary population.

They studied a metapopulation model with full dynamic in each patch in order to describe sexually transmitted diseases. They computed the basic reproduction number R 0 which is the average number of new infective generated by a single infective living in a completely susceptible population during his/her entire life span [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Diekmann | Mathematical epidemiology of infectious diseases. Model building, analysis and interpretation[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF][START_REF] Artzrouni | Transmission Probabilities and Reproduction Numbers for Sexually Transmitted Infections with Variable Infectivity: Application to the Spread of HIV Between Low-and High-Activity Populations[END_REF]. They also established bounds on R 0 . We analyze the global stability of the equilibria of this model. We use a theorem of [START_REF] Hirsch | The dynamical systems approach to differential equations[END_REF] and confirm the conjecture of [START_REF] Arino | A multi-city epidemic model[END_REF]. [START_REF] Arino | A multi-species epidemic model with spatial dynamics[END_REF] studied a multi-species multi-patch SEIR model, using the concept of asymptotic autonomous systems [START_REF] Castillo-Chavez | Asymptotically autonomous epidemic models[END_REF]) to prove the global stability of the disease free equilibrium (DFE). When R 0 > 1, they established the uniqueness and stability of the endemic equilibrium numerically. [START_REF] Arino | Quarantine in a multi-species epidemic model with spatial dynamics[END_REF] introduced quarantine in a SEIRS compartmental model for a multi-species multi-patch disease. [START_REF] Wang | An epidemic model in a patchy environment[END_REF] in a SIS model with "pseudo mass-action", assuming that the graph of migration is strongly connected, proved that the DFE is globally attractive when R 0 < 1, and that the disease is uniformly persistent when R 0 > 1. Jin and [START_REF] Wang | The effect of population dispersal on the spread of a disease[END_REF] showed that if R 0 > 1, there exists a unique endemic equilibrium which is globally asymptotically stable when the dispersal rates of susceptible and infective are equal. [START_REF] Auger | The Ross-Macdonald model in a patchy environment[END_REF] considered the Ross-Macdonald malaria model on n patches.

They proved that if R 0 ≤ 1, the DFE is globally asymptotically stable (GAS), and when R 0 > 1, the unique endemic equilibrium is also GAS. [START_REF] Iggidr | Analysis of a patchy sis epidemic model[END_REF] described the spatial propagation of a disease which confers no immunity. They proved that if R 0 ≤ 1, then the DFE is GAS; if R 0 > 1 there exists a unique GAS endemic equilibrium. [START_REF] Norman | A Host-Host-Pathogen Model with Vaccination and its Application to Target and Reservoir Hosts[END_REF] considered a multi-groups model in order to study the effectiveness of different vaccination strategies.

The migration model

Notations and definitions:

Residents of patch i: individuals who were born and reside in patch i;

Travelers: individuals who, at the current time, are not in the patch they reside in; n: the total number of patches;

N ij : the total number of residents of patch i who are present in patch j at time t; N r i = n j=1 N ij : the total resident population of patch i at time t; N p i = n j=1 N ji : the population of patch i at time t including both residents and travelers; g i > 0 : the per head rate at which residents of patch i leave it per time unit; m ij ≥ 0 : the proportion of individuals leaving patch i to patch j; g i m ji : the transfer rate from patch i to patch j; r ij ≥ 0 : the rate at which residents of patch i return to patch j; d : the death rate, equal to the birth rate.

We use the convention that m ii = 0 and r ii = 0, so that when g i > 0, n j=1 m ji = 1. Individuals do not give birth out of their resident patch and death occurs everywhere.

The migration model is the same as the one considered by [START_REF] Arino | A multi-city epidemic model[END_REF]:

N ii (t) = d (N r i (t) -N ii (t)) + n j=1 r ij N ij (t) -g i N ii (t), (1)
and for

j = i N ij (t) = g i m ji N ii (t) -r ij N ij (t) -d N ij (t). ( 2
)
We assume that some travelers return home. The return matrix R is R(i, j) = r ij and the outgoing matrix M is M (i, j) = g j m ji . We assume that these two matrices have the same zeros and nonzeros pattern, because they represent the return to i from j and the outgoing travel from j to i respectively. We assume that these matrices are irreducible, so that the n patches cannot be separated in two isolated groups. It is always possible to reduce the global study to that of irreducible components, thus our assumption does not reduce the generality of our results.

N r i (t) = 0 means that the population of residents in each patch is constant. However the population size in a patch i can vary.

Theorem 1 (Arino and van den Driessche ( 2003)) The system given by Eq. ( 1) and ( 2), for the initial value N ij (0) > 0, has the global asymptotically stable equilibrium

N ii = 1 1 + g i C i N r i (3)
N ij = g i m ji d + r ij 1 1 + g i C i N r i ( 4
)
where

C i = n k=1 m ki d + r ik .
The state of the system is given by n 2 components N ij ; it is an element of the nonnegative orthant [0, +∞) n 2 . Define the line-matrix N by: N = (N 11 , N 12 , . . . , N 1n , N 21 , . . . , N 2n , . . . . . . , N n1 , . . . , N nn ) T .

(5)

The migration model becomes

N (t) = M N (t), (6)
where M = diag(M ii ) is a block-diagonal matrix, with the block M ii for all i given by

M ii =             -g i r i2 + d r i3 + d . . . r in + d g i m 2i -r i2 -d 0 . . . 0
g i m 3i 0 -r i3 -d . . . 0 . . . . . . . . . . . . . . . g i m ni 0 . . . 0 -r in -d             . ( 7
)
3 The complete model

We assume that the transitions are independent of the disease status, which implies that infectious and susceptible individuals have the same migration rate. With this assumption the demographic change is given by Eq. ( 6). Let S ij (t) and I ij (t) be respectively the total number of susceptible and infective individuals in patch i who are present in patch j at time t: N ij (t) = S ij (t)+I ij (t) for all i, j = 1, . . . , n. We denote by β ijk > 0, the proportion of infectious contacts in patch j between a susceptible of patch i and an infective from patch k, and κ j > 0 the average total number of such contacts in patch j per time unit [START_REF] De Jong | How does transmission of infection depend on population size? In Epidemic models. Their structure and relation to data[END_REF][START_REF] Hethcote | The mathematics of infectious diseases[END_REF]. Subsequently, the infection rate of individuals in patch j is:

n k=1 κ j β ikj S ij (t) I kj (t) N p j (t) . (8)
The recovery rate γ is the same for all patches. The complete epidemic model is then:

I ii (t) = n k=1 r ik I ik (t) -g i I ii (t) + n k=1 κ i β iki S ii (t) I ki (t) N p i (t) -(γ + d) I ii (t) (9) S ii (t) = n k=1 r ik S ik (t) -g i S ii (t) - n k=1 κ i β iki S ii (t) I ki (t) N p i (t) + d (N r i -S ii (t)) + γ I ii (t) (10)
and for j = i,

I ij (t) = g i m ji I ii (t) -r ij I ij (t) + n k=1 κ j β ikj S ij (t) I kj (t) N p j (t) -(γ + d) I ij (t) (11) S ij (t) = g i m ji S ii (t) -r ij S ij (t) - n k=1 κ j β ikj S ij (t) I kj (t) N p j (t) -d S ij (t) + γ I ij (t). ( 12
) As S ij (t) = N ij (t) -I ij (t)
, it is sufficient to study Eq. ( 1), ( 2), ( 9), and ( 11). The differential equations governing the total number of infective individuals are:

I ii (t) = n k=1 r ik I ik (t) -g i I ii (t) + n k=1 κ i β iki S ii (t) I ki (t) N p i (t) -(γ + d) I ii (t) (13)
and for j = i,

I ij (t) = g i m ji I ii (t) -r ij I ij (t) + n k=1 κ j β ikj S ij (t) I kj (t) N p j (t) -(γ + d) I ij (t). ( 14
)
Define N p the vector of R n 2 given by:

N p = (N p 1 , N p 2 , • • • , N p n , N p 1 , • • • , N p n , • • • • • • , N p n ) T . ( 15
)
N p 0 means thats all the components of vector N p are positive. Using the same ordering as that of N ij , we define the vector I. The system given by Eq. ( 13) and ( 14) is written as:

I (t) = D I(t) -(γ + d) I(t) + diag(N p (t)) -1 diag(N (t) -I(t)) B I(t). ( 16
)
the block diagonal matrices D ii are defined by:

D ii (i, k) = r ik , D ii (k, i) = g i m ik , D ii (k, i) = -g i , D ii (k, k) = -r ik . (17)
namely:

D ii =             -g i r i2 r 13 • • • r in g i m 2i -r i2 0 • • • 0 g i m i3 0 -r i3 • • • 0 . . . . . . . . . . . . . . . g i m ni 0 • • • 0 -r in             . ( 18
)
Denoting e ij the canonical basis of R n 2 , the matrix B is defined by:

B e ij = n k=1 κ j β ikj e kj . (19)
B has the form

B =       B 11 • • • B 1n . . . . . . . . . B n1 • • • B nn       . ( 20
)
The matrix B is the block matrix B = (B jk ), where each block B jk is an (n × n) diagonal matrix. The (i, i) entry of B kk for every i is equal to κ i β kki ; for j = k the (i, i) entry of

B jk for every i is κ i β kji . For every k, B kk =       κ 1 β kk1 • • • 0 . . . . . . . . . 0 • • • κ n β kkn       , (21)
and for every (j, k), with j = k,

B jk =       κ 1 β kj1 • • • 0 . . . . . . . . . 0 • • • κ n β kjn       . ( 22
)
The complete SIS system is given by:

           N (t) = M N (t)
I (t) = D I(t) -(γ + d) I(t) + diag(N p (t)) -1 diag(N (t) -I(t)) B I(t). ( 23
)

Properties

We denote by "≤ "the point-wise ordering in R n , which is the ordering generated by the cone R n + . We also define the classical ordering in R n + by:

x < y if for any i x i ≤ y i and x = y x y if for any i x i < y i .

Positively Invariant Set

We define the vector:

N r = (N r 1 , N r 1 , • • • , N r 1 , N r 2 , • • • , N r 2 , • • • , • • • , N r n , N r n , • • • , N r n ) T . ( 24
)
Proposition 1 Define the set

K = {(N, I) : 0 ≤ N ≤ N r ; 0 ≤ I ≤ N r }.
Then K is a compact positively invariant set of Eq. ( 23).

Proof

We show that the vector field points inward from the faces of K.

If N ii = 0, then

N ii (t) = d N r i + n j=1 r ij N ij (t) ≥ 0. ( 25
) If ∀j = i, N ij = 0, then N ij (t) = g i m ji N ii (t) ≥ 0. If N ii = N r i , then ∀j, N ij (t) = 0 and N ii (t) = -d N ii (t) -g i N ii (t) ≤ 0. If ∀j = i, N ij = N r i , then N ii = 0 and N ij (t) = -r ij N ij (t) -d N ij (t) ≤ 0.
If I ii = 0, ∀i, then

I ii (t) = n k=1 r ik I ik (t) + n k=1 κ i β iki N ii (t) I ki (t) N p i (t) ≥ 0. ( 26
)
If ∀j = i, I ij = 0, then

I ij (t) = g i m ji I ii (t) + n k=1 κ j β ikj N ij (t) I kj (t) N p j (t) ≥ 0. ( 27
)
If

I ii = N r i , then I ij (t) = 0 ∀j = i, and I ii (t) = -g i I ii (t) -(γ + d) I ii (t) ≤ 0, because r ii = 0. If ∀j = i, I ij = N r i , then I ii (t) = 0, N ij -I ij (t) = N ij (t) -N r ≤ 0, and
I ij (t) = -r ij I ij (t) + n k=1 κ j β ikj (N ij (t) -N r ) I kj (t) N p j (t) -(γ + d) I ij (t) ≤ 0. ( 28
)
Hence the vector field associated with Eq. ( 1), ( 2), (13), and ( 14) points inward from the faces of K, so K is positively invariant under the flow of Eq. ( 23). Moreover, K is a closed bounded subset of R 2n 2 .

4.2 Reduction of Eq. ( 23)

Eq. ( 23) is a triangular system. Theorem 2 reduces the stability analysis to that of a smaller system.

Theorem 2 [START_REF] Vidyasagar | Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability[END_REF] Theorem 3.1 and 3.4) Consider the following

C 1 system      x (t) = f (x(t)) x ∈ R n , y ∈ R m y (t) = g(x(t), y(t)) (29)
with an equilibrium point (x * , y * ): f (x * ) = 0 and g(x * , y * ) = 0.

If x * is globally asymptotically stable (GAS) in R n for the system x (t) = f (x(t)), and if y * is GAS in R m for the system y (t) = g(x * , y(t)), then (x * , y * ) is locally asymptotically stable for Eq. ( 29). Moreover, if all trajectories of Eq. ( 29) are forward bounded, then (x * , y * ) is GAS for Eq. ( 29).

We consider Eq. ( 23) in the positively invariant compact set K and we know that the subsystem N (t) = M N (t) has a unique globally asymptotically stable equilibrium defined in Theorem 1. We denote this equilibrium by N , and the corresponding vector by N p . From Theorem 2 it is sufficient to study the stability of the reduced system: (2002). We define

I (t) = D I(t) -(γ + d) I(t) + diag(N p (t)) -1 diag(N (t) -I(t)) B I(t). (30
F = diag(N p ) -1 diag(N -I) B I, (31)

Result

Theorem 3 We consider Eq. ( 23) on K.

• If R 0 ≤ 1 then the system is globally asymptotically stable at the disease free equilibrium.

• If R 0 > 1 then there exists a unique endemic equilibrium (N , Ī) with I 0 which is globally asymptotically stable on

K \ {(N, 0), N ∈ [0, +∞) n 2 }.

Proof

It is sufficient to study the reduced Eq. ( 30) which is rewritten as:

I (t) = (F + V ) I(t) -diag(N p (t)) -1 diag(I(t)) B I(t). (36)
To prove this theorem, we use [START_REF] Hirsch | The dynamical systems approach to differential equations[END_REF] theorem:

Theorem 4 (Hirsch, 1984: Theorem 6.1) Let F be a C 1 vector field in R q , whose flow φ preserves R q + for t ≥ 0 and is strongly monotone in R q + . Assume that the origin is an equilibrium and that all trajectories in R q + are bounded. If the matrix-valued map DF : R q → R q × R q is strictly decreasing, in the sense that if x < y then DF (x) > DF (y) (37

)
then either all trajectories in R q + \ {0} tend to the origin, or there is a unique equilibrium p, (p 0) in the interior of R q + and all trajectories in R q + \ {0} tend to p.

The conclusion and the proof of this theorem are similar when R q + is replaced by a positively invariant subset

K of R q + .
Let

X(I) = (F + V ) I -diag(N p ) -1 diag(I) B I, (38)
the C 1 vector field in K.

The flow preserves K for t ≥ 0. The derivative DX is:

DX(I) = D -(γ + d) I n 2 + diag(N p ) -1 diag(N -I) B -diag(N p ) -1 diag(B I).
(39) DX(I) being an irreducible n 2 × n 2 Metzler matrix, the flow of X is strongly monotone in K. Because each row of B is nonnegative and nonzero, and the matrix-valued map DX is a decreasing function of I, so DX(I) strictly decreasing: if I 1 < I 2 then DX(I 1 ) > DX(I 2 ). Applying the Theorem 6.1 of [START_REF] Hirsch | The dynamical systems approach to differential equations[END_REF], we deduce that either all trajectories in K tend to the origin (which is the disease free equilibrium), or there is a unique equilibrium in the interior of K and all trajectories in K \ [0, +∞) n 2 × {0}

tend to this equilibrium.

The stability modulus α(M ) of a matrix M is the largest real part of the elements of the spectrum Spec(M ) of M .

α(M ) = max λ∈Spec(M ) Re(λ). ( 40
)
The Jacobian J(0) of Eq. ( 30) at the disease free equilibrium is J(0) = F + V . F ≥ 0 and V being a nonsingular Metzler matrix, F + V is a regular splitting of J(0). From [START_REF] Varga | Matrix iterative analysis[END_REF], ρ(-F V -1 ) < 1 is equivalent to α(F + V ) < 0. Hence the disease free equilibrium is locally asymptotically stable. Thanks to Hirsch's theorem, the origin is then globally

asymptotically stable if R 0 < 1. If R 0 > 1, this is equivalent to α(J(0)) > 0.
The disease free equilibrium is then unstable; there exists a unique attracting endemic equilibrium I 0. This endemic equilibrium I satisfies:

(D -(γ + d)) I + diag(N p ) -1 diag(N -I) B I = 0. ( 41
)
Thanks to the non-negativeness of B and to the fact that I 0 we get:

DX(I) I = -diag(N p ) -1 diag(B I) I < 0. ( 42
)
Using the fact that DX(I) is a Metzler matrix, Eq. ( 42) implies that it is stable (Berman and Plemmons, 1994: criterion 28 of Theorem 6.2.3).

The stability modulus then satisfies α(DX(I)) < 0. This proves the local asymptotic stability of I, and using Theorem 4 we deduce that I is globally asymptotically stable if

R 0 > 1.
To complete the proof, we consider the case R 0 = 1, which is equivalent to α(F + V ) = 0. As F + V is an irreducible Metzler matrix, there exists a positive vector v such that (F + V ) T v = 0. We consider the following Lyapunov function

V (I) = I| v , (43)
on K. The derivative along the trajectories is:

V = -diag(N p ) -1 diag(I) BI | v ≤ 0. ( 44
)
This proves the local stability of the disease free equilibrium. By the Theorem 4 of Hirsch we are necessarily in the case where the disease free equilibrium is attractive, which means that the disease free equilibrium is globally asymptotically stable.

Simulations

We consider two patches with parameters chosen so as to describe gonorrhea [START_REF] Arino | A multi-city epidemic model[END_REF]. The time unit is a day; the recovery rate is γ = 0.04; the average lifespan is 75 years, and d = 1/(75 × 365). We assume g 1 , g 2 > 0 so that m 12 = m 21 = 1. The disease transmission coefficients are equal in each patch for all contact types: β ikj = β k , ∀i, j = 2, and

F = 1 N p 1 N p 2          κ 1 β 1 N 11 N p 2 0 κ 1 β 1 N 11 N p 2 0 0 κ 2 β 2 N 12 N p 1 0 κ 2 β 2 N 12 N p 1 κ 1 β 1 N 21 N p 2 0 κ 1 β 1 N 21 N p 2 0 0 κ 2 β 2 N 22 N p 1 0 κ 2 β 2 N 22 N p 1          = 1 N p 1 N p 2    F11 F11 F22 F22    . ( 45
) -V =          g 1 + γ + d -r 12 0 0 -g 1 r 12 + γ + d 0 0 0 0 r 21 + γ + d -g 2 0 0 -r 21 g 2 + γ + d          =    V 11 0 0 V 22    . (46)
V

-1 11 = 1 ∆ 1    r 12 + γ + d r 12 g 1 g 1 + γ + d    = 1 ∆ 1 V -1 11 . (47)
V

-1 22 = 1 ∆ 2    g 2 + γ + d g 2 r 21 r 21 + γ + d    = 1 ∆ 2 V -1 22 . ( 48
) -F V -1 = 1 ∆ 1 ∆ 2 N p 1 N p 2       ∆ 2 F11 V -1 11 ∆ 1 F11 V -1 22 ∆ 2 F22 V -1 11 ∆ 1 F22 V -1 22       , (49)
where ∆ 1 = (γ + d) (γ + d + g 1 + r 12 ) and ∆ 2 = (γ + d) (γ + d + g 2 + r 21 ).

The matrix of Eq. ( 6) is a 4 × 4 matrix with rank two. Proposition 3 will help us to reduce the computation of its spectral radius to that of a 2 × 2 non singular matrix.

I n 0 CA -1 I n       A 0 0 D -CA -1 B       I n A -1 B 0 I n    .
Let λ ∈ Spec(P ), then The application of Proposition 3 to (-F V -1 ) yields:

R 0 = 1 ∆ 1 ∆ 2 N p 1 N p 2 ρ F11 ∆ 2 V -1 11 + ∆ 1 V -1 22 F22 F -1 11 . Or R 0 = 1 ∆ 1 ∆ 2 N p 1 N p 2 ρ(∆ 2 V -1 11 F11 + ∆ 1 V -1 22 F22 ). (54)
With Eq. ( 54), the calculation of R 0 is straightforward because we deal with a 2 × 2 non singular matrix.

Set

A = ∆ 2 V -1 11 F11 + ∆ 1 V -1 22 F22 . Then R 0 = tr(A) + (tr(A) 2 -4 det(A)) 1/2 2 ∆ 1 ∆ 2 N p 1 N p 2 . ( 55
)
We simulate in order to underline the effect of the migration of individuals on the metapopulation. We compute R 0 using Eq. ( 55). The first parameters are :

κ 1 = 1, κ 2 = 1, 1 = g 2 = 0.05, β 1 = 0.02, β 2 = 0.05, r 12 = 0.2, r 21 = 0.05, N r 1 = 1500, N r 2 = 1500, γ = 1/25, d = 1/(75 × 365).
When the two patches are isolated, the disease will die out in patch 1 with R 1 0 = 0.5 while it is endemic in patch 2 with R 2 0 = 1.25. When the two patches exchange people according to the parameter values, R 0 = 0.96: the solutions go to zero and the disease dies out in the metapopulation (Figure 1(a)).

In Figure 1

) 4. 3

 3 Basic reproduction ratio R 0The basic reproduction ratio was computed by[START_REF] Arino | A multi-city epidemic model[END_REF]. Here, we express it in terms of the matrices B and D. We use the framework of[START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations[END_REF],[START_REF] Diekmann | Mathematical epidemiology of infectious diseases. Model building, analysis and interpretation[END_REF], and van den Driessche and Watmough

  Schur complement and the properties of the determinant. Recall that, if A is an invertible matrix, and M = bloc matrix, then the Schur of A in M is the matrix D -CA -1 B. The Schur complement is usually used to obtain the following factorization: M =   

  (b), we keep the same parameters as in Figure 1(a) except that r 12 is decreased significantly from 0.2 to 0.02. The reproduction number in the metapopulation grows to R 0 = 1.1, whereas R 1 0 and R 2 0 remain unchanged because they are not affected by a modification of the mobility coefficients g i , m ij , and r ij . The disease becomes endemic in the metapopulation. Decreasing r 12 accounts for lengthening the stay of individuals in patch 2 where the disease is initially endemic. The people from patch 1 spend most of their time in the endemic patch 2. The same explanations hold when in Figure 1(a) g 1 increases from 0.35 to 0.85 (R 0 = 1.03) because g 1 reflects the flow of individuals from patch 1 to patch 2. This is shown in Figure 2(c).

Figure 3

 3 Figure 3(a) shows that the metapopulation remains disease free when all patches are initially disease free. In Figure 3(a), we keep the parameters of Figure 1(a) and decrease β 2 from 0.05 to 0.03.

Figure 3

 3 Figure 3(b) shows that the metapopulation remains endemic when all patches are initially endemic. This is the consequence of our Theorem 3 and of the Theorem 6 of Arino and van den Driessche (2003). In Figure 3(b), we keep the parameters of Figure 1(a)

Figure 2 :

 2 Figure 2: (c) The solutions of the system (23) converge to the endemic equilibrium when the parameters are the same as in Figure 1(a) except g 1 = 0.85 instead of g 1 = 0.35 : R 1 0 < 1, R 2 0 > 1, R 0 > 1.

Figure 3 :

 3 Figure 3: (a) The solutions of the SIS system (23) converge to DFE when all patches are initially disease free. The parameters are g 1 = 0.35, g 2 = 0.05, β 1 = 0.02, β 2 = 0.03, r 12 = 0.2 : R 10 < 1, R 2 0 < 1, R 0 < 1. (b)The solutions of the SIS system (23) converge to the endemic equilibrium when all patches are initially endemic. The parameters are g 1 = 0.35, g 2 = 0.05, β 1 = 0.08, β 2 = 0.05, r 12 = 0.2 : R 1 0 > 1, R 2 0 > 1, R 0 > 1. 24

  Moreover, if P 12 commutes with P 22 , ρ(P ) = ρ (P 11 + P 22 ) .

				
	Proposition 3 Let P =	 	P 11 P 12 P 21 P 22	  be a rank n matrix with P 11 , P 12 , P , P 22 n × n
	square matrices.
	If P 12 is invertible, and P 12 P 21 -P 12 P 22 P -1 12 P 11 = 0, then
			ρ(P ) = ρ P 11 + P 12 P 22 P -1 12 .	(50)

Acknowledgments

We thank the anonymous referees for their valuable comments and suggestions.

[Figure 1 about here.] [Figure 2 about here.] [Figure 3 about here.]

 [START_REF] Arino | A multi-city epidemic model[END_REF]showed local stability and studied global behavior numerically. They conjectured the global stability of equilibria depending on the basic reproduction number. We turned this into Theorem

the function of appearance of new infections in infectious compartments and

the transfer in compartments by all other means.

The Jacobian F of F at the disease free equilibrium is:

and the Jacobian V of V is:

The matrix D is a Metzler matrix (nonnegative off diagonal terms) with a zero column sum. This implies that 0 is a simple eigenvalue of D, the other eigenvalues having a negative real part [START_REF] Jacquez | Qualitative theory of compartmental systems[END_REF].

The eigenvalues of V are the ones of D subtracted of (γ + d), hence V is a stable Metzler matrix with -(γ + d) as stability modulus. This implies that V is nonsingular.

We then apply the results of van den [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] to obtain the basic reproduction ratio:

Proposition 2 The basic reproduction ratio is given by

where ρ(A) denotes the spectral radius of the matrix A.