Mathematically modelling competitive ion absorption in a polymer matrix

A recent publication [Lindén et al., RSC Adv., 2014, 4, 25063–25066] describing a material with extreme efficiency and selectivity for adsorbing copper opens up the possibility of large scale purification technologies and new materials for prevention of biological growth, both in bulk and on surfaces. On the molecular scale, the material’s efficiency and selectivity depends on competitive metal ion diffusion into a nanometer thin polymer film and on competitive binding of metal ions to specific ligand sites. We present and explore two reaction-diffusion models describing the detailed transport and competitive absorption of metal ions in the polymer matrix studied in [Lindén et al., RSC Adv., 2014, 4, 25063– 25066]. The diffusive transport of the ions into an interactive, porous polymer layer of finite thickness, supported by an impermeable substrate and in contact with an infinite electrolyte reservoir, is governed by forced diffusion equations, while metal ion binding is modelled by a set of coupled reaction equations. The qualitative behavior observed in experiments can be reproduced and explained in terms of a combination of (a) differing ion diffusive rates and, effectively, (b) an ion exchange process that is superimposed on independent binding processes. The latter’s origin is due to conformational changes taking place in the polymer structure in response to the binding of one of the species. The results of simulations under a diverse set of parameter conditions are discussed in relation to experimental observations.


Introduction
Water contaminated with toxic metal ions is oen a result of industrial processes 2 but also occurs through antifouling paints leaching metal-based biocides. 3,4Increased awareness of the toxic effects and economic aspects of reuse of water, sludge and sediment and recovery of metal waste are driving forces for research in new and improved metal-removal technologies 2,5 .A wide range of techniques has been studied over past decades, including chemical precipitation, ion-exchange, adsorption, membrane ltration, coagulation-occulation, otation and electrochemical methods. 6dsorption and ion-exchange technologies have proven to be effective methods for the removal of various metals from aqueous solutions. 6,7][8] The separation process of the metal ions from the aqueous phase is typically described as adsorption onto the surface of the solid phase of either the adsorbent 9 or the ion-exchange resin. 2 Wan Ngah et al. 9 evaluated the adsorption capacities and rates of Cu(II) ions onto chitosan and chitosan blended with poly(vinyl alcohol) (PVA).They found that the Langmuir model correlated well with the adsorption isothermal data and that the kinetic experimental data correlated well with the second-order kinetic model which would indicate the rate-limiting step as chemical sorption.Beatty et al. 2 also analysed the performance of their ion-exchange resin using a Langmuir isotherm and second-order kinetic model.They described their slower rate of adsorption in comparison to a commercially available resin as being dependant on a longer ion diffusion path into the polymer layer on the silica support material.
Although adsorption isotherms such as the Langmuir isotherm play an important role in the search for an ideal metal ion adsorbent, they might face limitations when the metal ion uptake is not solely adsorption. 10In one study of chitosan and its glutaraldehyde cross-linked derivative it was concluded that the metal ion uptake mechanism was absorption within the polymers rather than adsorption on the polymer surface. 11This indicates the importance of investigating the possibility of considering the metal ion uptake of a pure polymer sorbent or polymer modied sorbent as an absorption process.Furthermore, the effect of conformational changes within polymers could also play an important role in the metal uptake mechanism. 12,13e address the following problem: to model the temporal and spatial development of competitive absorption of metal ions within a polymer matrix.The example system we have in mind to model is a polyethyleneimine (PEI) layer adsorbed onto an impermeable planar surface and exposed to a bulk electrolyte mixture.Recently, Lindén and coworkers 1 experimented with a surface-bound, cross-linked PEI layer exposed to either an articial sea water solution or actual sea water.They found that the PEI layer exhibited preferential absorption for certain metal cations, namely copper (Cu 2+ ), zinc (Zn 2+ ) and cadmium (Cd 2+ ), in that order of preference.The competitive absorption of metal ions into the PEI layer was performed at an equimolar concentration of about 3 mM, 1 which is roughly 200 times higher than the levels at which the ions are present in sea water. 14,15][18] Linden 1 has recently shown that polyethyleneimine (PEI) is extremely effective and selective to copper in salt and fresh water and that it self-adsorbs into a thin lm on porous silica particles.The latter allows for very large surface area exposure.Furthermore, the cross-linked polymer is mechanically and chemically stable in a wide pH range.Taken together, the material properties of PEI and porous silica particles allow for large scale technologies for purifying water from copper as well as preventing biological growth at surfaces and in bulk.Given the need for copper purication and the large scale potential of the PEI/silica particle system we aim to increase our understanding of the molecular-scale processes by addressing the following problem in more detail: to model the temporal and spatial development of competitive absorption and chemical binding of metal ions within a polymer matrix.
0][21] Rivas et al., 20 working with PEI tetramers (TETA) in solution, showed that the TETA molecules could take up copper, and Kobayashi with coworkers in ref. 19, who also worked in solution, found that the branched structure of PEI gave a more favourable complexation than the linear structure.7][18] However, what was not explored prior to the work of Lindén et al. 1 was the competitive absorption behavior of PEI in the presence of a solution mixture, the predominant results of which are shown in Fig. 1.In particular, Lindén et al. 1 found that the cross-linked PEI matrix displayed an anomalous preference for copper in the following sense.At a very early stage post exposure, zinc ions, and to a lesser extent cadmium ions, were rapidly absorbed forming complexes with specic sites in the matrix.Later, increasing amounts of copper complexes were detected.However, rather than each metal ion complex increasing in number monotonically with prolonged exposure to a bulk electrolyte reservoir, it was found that the copper ion complexes increased signicantly, reaching a steady state plateau value, while the zinc and cadmium complexes reached a maximum amount and then decreased to levels below instrument sensitivity.Effectively, the quantities of bound zinc and cadmium approached zero at steady state.It is this entire time series of phenomena that we seek to address in this paper.Reasons for the ndings were postulated, as were the possible applications of the phenomena, such as the extraction of minerals from seawater (specically, the extraction of copper from the oceans), the renement and purication of industry waste water, as well as the more community-based treatment of recycled water. 1 In explanation for the two phases of uptake, Lindén et al. 1 stated that "the binding kinetics were faster for zinc but with copper binding being thermodynamically favored".On the face of it, this would seem a natural conclusion to draw.However, the physical system proves to be more complex.In this paper we present and investigate quantitatively two alternate scenarios which explain the competitive behavior, at least qualitatively.In the rst scenario we postulate that the differential binding phenomenon is simply not due to faster binding kinetics on the part of Zn 2+ and a steady state or thermodynamic preference toward Cu 2+ binding, but rather to a difference in the diffusive properties of the two ions in conjunction with an ion exchange reaction.In the second scenario, in recognition of the ndings of Nurchi et al., 13 who argued that the binding of Cu 2+ induces a conformational change in the matrix structure, we postulate that the competitive binding phenomenon is due to a conjunction of differential diffusion and a conformational change upon Cu 2+ binding that both precludes the binding of Zn 2+ ions as well as ejecting already bound Zn 2+ in the neighborhood of the copper binding site.
Although the (two) models we introduce and explore below target the aforementioned experimental system directly, it is not unreasonable to anticipate that the models and the simulation results derived from them have applicability in other contexts, notably in biology, dealing, for instance, in ion specic membrane transport. 22he remainder of this paper is organized as follows.In Section 2 we describe the two scenarios hypothesized to explain the ndings, as well as introducing other essential particulars of the models.In Section 3, the mathematical descriptions of the physical models are outlined.Here, the key model parameters are introduced.In that section we also outline the numerical solution method followed.Section 4 presents some illustrative example simulations, which explore the important regions of parameter space.We summarize our ndings and provide concluding remarks in Section 5.

Polymer matrix geometry & kinetic transport mechanisms
In line with the experimental system, we adopt the model of a thin, charged polymer matrix of thickness D that has been adsorbed to an uncharged, solid and impermeable substrate.We assume in the following that D remains a constant and is not responsive to ion uptake (an alternate model approach showing this effect can be found in ref. 23 and 24).We further assume the polymer matrix to be homogeneous and uniformly porous.Moreover, it is homogeneously permeable to mobile ions.We assume the matrix to possess a constant (total) density of binding sites, n PEI tot , which are either equally accessible to both binding cations (Scenario I) or are equally divided into a number that are accessible to only Cu 2+ (n PEI 1,tot ) or only Zn 2+ (n PEI 2,tot ) (Scenario II).In either case, the binding sites are assumed to be negatively charged with unit valence (z PEI ¼ z PEI 1 ¼ z PEI 2 ¼ À1).The adsorbed layer is in contact with a 2 : 1 mixed bulk electrolyte comprising Cu 2+ , Zn 2+ and Cl À ions.
The important processes taking place in this problem include mobile ion diffusion into or out from the polymer matrix, with ions sourced from the electrical double layer adjacent to the polymer matrix.In addition, there are two fundamental chemical binding reactions, one for each of the cations Cu 2+ and Zn 2+ and their respective binding sites within the matrix.The details of these intrinsic reactions are given in the next section.In Section 4 we show that these effects alone are insufficient to explain the experimental observations.Therefore, in addition to the core binding reactions we invoke a secondary reaction mechanism, which will be explained in Section 3 via two model formulations, that completes the picture and quite convincingly accounts for observed phenomenon.Finally, the outcome of all the binding events (as well as the underlying charged state of the matrix) results in an inhomogeneous distribution of charge that establishes an inhomogeneous electric potential eld.The eld gradient in turn inuences the rate of inux of ions (indeed it inuences which ions enter the matrix as well as their diffusion rates).
With respect to the electric eld, we envisage that the mobile ions external to the matrix diffuse sufficiently fast that the ionic double layer outside the matrix is always in a state of quasi equilibrium.That is, the ions are sufficiently responsive on the time scale of the processes taking place within the layer that the double layer adjusts immediately in response to changes in the state of polymer charge.At the low values of accumulated charge and low bulk electrolyte concentrations considered, we assume that the electrical double layer outside the matrix is adequately described by the Poisson-Boltzmann model. 25,26We remark that although we will associate differences in ion diffusivities with differences in effective ionic radii (see Section 4), our models do not include any explicit ion size effects, which could not only inuence the double layer structure adjacent to the polymer layer at high concentrations, 27 but would likely play an inuential role in shiing the balance in the competition for absorption sites even at low electrolyte concentrations.We envisage the latter effect to be akin to the site exclusion phenomenon simulated in the hydrated ion model of competitive ion adsorption. 28We reserve for the future the task of considering this effect, whose signicance is comparable to and indeed warrants a more accurate description of the polymer structure itself, which is presently unavailable.

Governing equations
The charge distribution within the polymer matrix comprises free ions, here denoted by the index i ¼ 1,.,n and bound charges, here denoted by the index j ¼ 1,.,m.As our focus is on the competitive binding of Cu 2+ and Zn 2+ in a polymer matrix, we restrict attention to determining the distribution of these two ions plus that of a common counterion, Cl À .Thus, n ¼ 3.In this same context, the bound species are PEI À , CuPEI + and ZnPEI + , as explained below.Thus, m ¼ 3 or 4, as explained below.The charge density is given as, rðx; tÞ ¼ X n i¼1 z i en i ðx; tÞ þ X m j¼1 z j en j ðx; tÞ; where n i and n j refer to the number densities of the respective species, i and j at position x at time t, with z i and z j being their respective valencies; e is the electronic charge.

Competitive binding scenario I: Ion exchange model
To describe the competitive absorption of Cu and Zn ions, we model the ion binding event here assuming a single type of binding site, denoted PEI À .In this case we rst identify two binding reactions, one for each of the ions.These reactions involve a mix of free ions and bound charge sites: Thus far, these reactions take place explicitly independent of each other and occur at different rates given by where we have denoted the molar concentration of component X by the symbol [X] (M); molar concentrations are related to number densities (cm À3 ) by the relation, n X ¼ N A 10 À3 [X], where N A is Avogadro's number, 6.022 Â 10 23 mol À1 .The rate constants of forward (+) and reverse (À) reactions are denoted by the parameters k AE .In addition to these binding reactions we include in this scenario the following exchange reaction, which occurs at the rate The rates, r 1 and r 2 , describe the rates of production of CuPEI + and ZnPEI + , respectively, as well as the consumption of Cu 2+ and Zn 2+ , the latter being the free mobile ions.Given the additional exchange reaction, the rate equations ( 4) and ( 5) must be supplemented by the effect of the exchange rate r 3 : and, similarly, These equations are explicit, ordinary differential equations in time, but are also implicitly functions of position x through the dependencies of the respective concentrations on position.That is, these equations must be solved locally (at x) throughout the polymer matrix.
This binding model must satisfy the conservation condition that the total number of sites remains xed irrespective of whether there is binding or not.The relevant equation is

Competitive binding scenario II: Polymer relaxation model
We formulate this model situation by assuming, at least initially, that there are two types of binding sites (labelled with subscripts 1 and 2), one for each of the divalent metal ions, Cu 2+ and Zn 2+ .In addition to the standard binding reaction equations, we presume here that site type 2 can exist in one of two binding states, an active binding state, labelled {2, a}, and an inactive binding state, labelled {2, ia}.The fundamental reaction equations are then slightly modied versions of eqn ( 2) and (3), where, apart from the numerical subscripts, there now appears in eqn (12) the reference {2, a} as a reminder that the chemical reaction takes place on an active type-2 binding site only.
A point of major differentiation between this scenario and that of Scenario I is the assumption that the binding of Cu 2+ to its specic binding site induces a local conformational change in the polymer matrix which converts an active type-2 binding site to its inactive binding state.This conversion, which may be structural or chemical, either prevents the binding of Zn 2+ to its specic binding site or induces an unbinding of an existing bound Zn 2+ .For the latter possibility, we imagine the two-step reaction process followed by a Zn 2+ unbinding or dissociation reaction Once in the inactive state, the forward reaction in eqn ( 14) is considered to be sufficiently rapid as to be a spontaneous event.Hence, one would expect k + 4 [ k À 4 .Indeed, assuming the limiting case of complete forward reaction, we can summarize eqn ( 13) and ( 14) in one reaction equation, The rate of production of free Zn 2+ (assuming complete dissociation) is then governed by the rate at which reaction (15)  proceeds, which is given by which is also the rate at which ZnPEI + 2,a is converted to PEI À 2,ia .To be completely consistent in our reasoning, we should actually not consider the possibility of a reverse reaction in eqn (16).That is, we should set k À 3 h 0 from the outset.Nevertheless, we include the reverse terms in eqn ( 15) and ( 16) for formal completeness only and point out that in our calculations we avoid adopting signicant values of k À 3 or k À 4 and so maintain the forward balance of the reactions.
As for the case of converting the state of a vacant type-2 binding site, we invoke the model reaction which occurs at the rate In all probability, the latter rate constants should satisfy the inequality k + 5 [ k À 5 , implying that the reverse reaction all but never takes place.Note also that based on our intuitive understanding of the conformational change process, we expect the k + 5 value to be different, but perhaps not signicantly different, from the k + 3 value.As in the case of the single site binding model, this two-site binding model must satisfy conservation conditions: the total number of sites remains xed irrespective of whether there is binding or not.The relevant equations are It is important to point out that Scenario I and II are different in two major respects.Firstly, Scenario I describes a process of physical exchange of Cu 2+ for Zn 2+ , with no alteration of the polymer matrix, while Scenario II describes a conformational reaction within the polymer resulting from the binding of Cu 2+ .These are quite distinct phenomena.The second point of differentiation is that Scenario I assumes a common type of binding site, equally accessible to both Zn 2+ and Cu 2+ , while Scenario II assumes specic sites for the respective ions.It is not difficult to imaging other scenarios that are intermediate or hybrids of the two cases.However, it is not the purpose of this paper to establish which of these is most applicable to the given experimental system.Our purpose is only to demonstrate that an additional interaction involving the ions (such as described by Scenarios I and II) is necessary to produce behavior reminiscent of the experiments.

Ion diffusion
The intrusion of ions into the polymer matrix is assumed to take place by means of (forced) diffusion.Exposure of the matrix to an electrolyte solution with ions at a nite bulk concentration and thus a larger chemical potential induces the transport of ions into the polymer layer.The process is largely driven by concentration gradients but also by an electric eld gradient arising from the inhomogeneous charge distribution.Three ions are assumed sufficiently mobile to diffuse into the matrix: Cu 2+ , Zn 2+ and Cl À .The rates at which ions diffuse into the matrix are here governed by their intrinsic diffusion constants and their respective ionic conductivities.However, it is well known that diffusion of dissolved molecules is reduced, relative to bulk diffusion, by polymer chains through several mechanisms. 29While it is necessary for a quantitative evaluation with experiments to consider these effects, for the present task of illustrating qualitative behavior, and in the absence of detailed quantitative experimental information about the present system, we shall assume, unless otherwise stated, that the values of the ion diffusion parameters are equal to their bulk values. 30The equations governing forced ion diffusion are where the terms on the le hand sides of the respective equations combined with the rst terms on the right hand sides conjointly comprise the usual free diffusion equations for the mobile ions, Cu 2+ , Zn 2+ and Cl À .The second terms on the right hand sides represent the leading order, nonideal activity, contributions to the gradient of the chemical potentials.In the mean-eld model these contributions involve the gradient of the self-consistent mean-eld electric potential.The nal terms on the right hand sides of the equations for Cu 2+ and Zn 2+ refer to the net consumption/production of these ions as determined by their respective sets of chemical absorption reactions.Absorption of the ions means they are no longer free to participate in the diffusive process.The effects, represented by the functions P Cu and P Zn , are therefore treated as that of sinks.The opposite is true of ions produced through the reaction events.The latter effects are treated as that of sources.In eqn (21), the l i parameters are ionic drag coefficients, while the D i parameters are ion diffusion constants.
Eqn (21) are complemented by boundary and initial conditions.Initially, we suppose there are no mobile ions present.Thus, At x ¼ 0 we suppose there is no ux of ions through the substrate while at x ¼ D we suppose there is a ready supply of mobile ions equal to the concentration of ions in the double layer adjacent to the polymer interface,

Electric potential
Given that the bound ions are also charged contributions, the concentrations will also need to be incorporated into the expression for the charge density.The explicit expression for the charge density is z j en j ðx; tÞ The non-uniform charge content of the polymer layer generates an electric eld which inuences both the rate and extent of ion intrusion into the polymer layer.The electric potential associated with that eld satises Poisson's equation, and, nally, the conditions of continuity of electric potential and electric displacement at x ¼ D, the interface between the polymer matrix and the electrolyte, 8 < : where, in the last equation, we have assumed that the matrix and bulk electrolyte have the same relative dielectric permittivity.
Although the equations for the electric potential depend explicitly only on the spatial variable, x, the potential itself is also implicitly dependent on time, t, through its dependence on the charge density which depends explicitly and implicitly on time.Consequently, Poisson's equation and the Poisson-Boltzmann equation need to be solved self-consistently throughout the time period of interest.Incidentally, the application of the Poisson-Boltzmann equation in the bulk solution implies the assumption of rapid equilibrium in the bulk.That is, irrespective of the time development within the polymer matrix, the electrolyte in the adjacent electrical double layer is taken to relax sufficiently fast to always be in a state of (quasi-) equilibrium given the prevailing conditions within the matrix.

Numerical solution approach
To solve the coupled systems of equations we rst nondimensionalise all variables and equations using the time and spatial scales of D m k 2 and k À1 , respectively, where We then introduce dimensionless dependent and independent variables and related parameters: where is the limiting ion conductance of species i, while l i is the ionic drag coefficient for that ion.
The equations to be solved are then written in dimensionless form.For example, non-dimensional versions of the diffusion equations take the form with a similar equation for Cl À , but without the last term on the right hand side.The Poisson equation becomes z l C l ðy; sÞ; 0\y\K: (34) Finally, the reaction equations (e.g., eqn (2) and ( 3)) typically appear in the form Note that some dimensionless reaction rate constants are not dened as per eqn (32), depending on their original dimensions.Nevertheless, dimensionless reaction rates are implicit in all of our working.We point out that for convenience we subsequently refrain from quoting the units of the different reaction rate constants.Instead, we document here that the units of k + 1 , k + 2 , k + 3 , k À 3 (Scenario I and II), k À 4 and k AE 5 are cm 3 mole À1 s À1 , while the units of k À 1 , k À 2 and k + 4 are s À1 .The units of k À 3 (rapid dissociation version of Scenario II) are cm 6 molec À2 s À1 .
It is possible in principle to consider all three sets of equations (dimensional or nondimensional versions) as subsystems of a larger system of partial differential equations (pdes) involving both time (t or s) and space (x or y) derivatives.Those equations not possessing explicit time (space) derivatives can be considered as pdes with zero time (space) derivative coefficients.The complete system could theoretically be solved as one selfconsistent system.Unfortunately, this turns out to be problematic in practice.In this work, we have adopted an iterative approach whereby we cycle through the solution of each subsystem in turn, assuming that the dependent variables not being solved for explicitly are inhomogeneous inputs in that calculation.The outputs of each subsystem calculation are then used as inputs for the solution of the other subsystems.The entire cycle of calculation is iterated until we achieve selfconsistent results for a given set of initial parameter conditions.The code developed was written in the Matlab® environment utilizing some of the inbuilt subroutines.The subsystem of reaction equations of the form in eqn (35) was solved using the function routine ode45, while the set of diffusion equations, eqn (33), was solved using the solver pdepe.With regard to eqn (34), this equation was integrated explicitly to give the result fðy; sÞ ¼ with F 0 to be determined from satisfaction of the continuity conditions at y ¼ K: In eqn (37), the right hand side of the second of the continuity equations is the derivative of the Poisson-Boltzmann solution: which explicitly involves the solution evaluated at y ¼ K. Eqn (38) and the rst of eqn (37) represent a single closed transcendental equation that can be solved, say iteratively, for F 0 .However, it turns out that for the charge conditions of interest here, the potential is sufficiently small in magnitude to justify invoking the linear Poisson-Boltzmann solution approximation outside of the polymer region.This greatly simplies the procedure allowing F 0 , and thus the entire potential function within the matrix, to be determined explicitly.In this case we nd that 4 Results and discussion

Scenario I
In this section we shall focus on Scenario I as an example of a system exhibiting the type of qualitative behavior found in the experiments of ref.
1.The working hypothesis of this scenario is that although both ions have a similar binding affinity for the PEI binding sites (or possibly that Cu 2+ has a greater affinity for binding than Zn 2+ ), the Zn 2+ metal ion has a more rapid diffusive entry into the polymer matrix and therefore has an earlier binding association with the polymer than does Cu 2+ .However, the hypothesis then proposed is that once in, the Cu 2+ ions preferentially replace already bound Zn 2+ ions to either a signicant fractional extent or completely, depending on the respective binding constants.Clearly, the important parameters to focus on in order to explore and test the hypothesis are the diffusion constants of the respective ions and the binding reaction rate constants introduced in the previous section.To this end (and so as to limit the very large parameter space inherent in the model) we focus attention on the principal parameters: k AE 1,2 , k AE 3 and D 1,2 , keeping other quantities constant.Quantities such as the thickness of the polymer matrix, w, the bulk concentrations of free electrolyte ions, C N Cu , C N Zn , C N Cl and the density of binding sites serve only to adjust the total number of bound Cu 2+ and Zn 2+ ions (i.e., the steady state) for a given set of our designated principal parameters.It is possible that under extreme parameter conditions, these secondary parameters may inuence even the qualitative behavior.However, it will become clear that under quite reasonable parameter conditions, the model described in Section 3 does support the hypothesis.Table 1 summarizes the xed parameters used in these calculations.
It has been shown 30 that the diffusion constants of the mobile electrolyte ions vary, depending on solution conditions.The ranges of diffusion constants, 0.5-1.0Â 10 À5 cm 2 s À1 and 2-5 Â 10 À6 cm 2 s À1 , determined experimentally for Cu 2+  We rst establish the reference case of ion binding, assuming all but equal diffusion constants and equal binding affinities.Under such conditions there is identical time and space development of the density of CuPEI + and ZnPEI + complexes within the polymer layer.One should note that steady state densities are obtained for all mobile ions well within the non-dimensional time frame adopted.The steady state concentrations of CuPEI + and ZnPEI + are dictated solely by the respective bulk ion concentrations and the underlying density of binding sites.The equivalence of the two ions results in negligible difference in either diffusive or absorptive behavior.This is summarized by the coincident curves in the centre of Fig. 2. The gure shows the time evolution of the integral amounts of CuPEI + and ZnPEI + in the matrix.Mobile metal ions occupying space in the matrix are not included in these integral quantities.The curves have converged to the theoretical limit depicted by the dashed line, which represents the case of equal share of binding sites to Cu and Zn.The separated curves (blue solid curves for CuPEI + and red broken curves for ZnPEI + ) on either side of the coincident pair depict the result of delaying the inux of Cu 2+ by reducing this ion's diffusion constant.The two values of R Cu ¼ 3.5 Å and 4.5 Å show progressively reduced Cu 2+ uptake with corresponding increased Zn 2+ uptake.The accumulating ZnPEI + , however, reaches a maximum and then decreases as CuPEI + continues to increase.However, the behavior is not due to any extraneous exchange process as negligible exchange has been assumed here (k AE 3 ¼ 1.0).Instead, what we observe is a natural approach of the two corresponding curves to their equal steady state values shown by the black dashed line in the gure centre.All things being equal, the differences in diffusive rates only affect the approach to equilibrium, not the equilibrium state itself.
Thus, although we have not been able to track numerically the time development further, we surmise that the two sets of curves will each ultimately asymptote to the theoretical nondimensional limit of 0.3153.(Numerical difficulties were experienced in solving the coupled system of diffusion equations, presumably because the faster species achieves steady state values in signicant regions of the matrix before the overall system reached steady state.We aim to address this problem mathematically in future work.) As our interest is with competitive absorption of metal ions over and above the intrinsic process depicted by the curves in Fig. 2 and governed by eqn (2) and (3), we show in Fig. 3 the time development of total metal ions bound in the matrix under equivalent diffusive and equal binding conditions, but with increasing inuence of the forward exchange reaction.Three cases are demonstrated: k + 3 ¼ 1.0, 1.0 Â 10 9 and 1.0 Â 10 12 .The blue solid curves again follow the accumulation of bound Cu 2+ throughout the polymer, while the red broken lines follow the accumulation of Zn 2+ .A nine order of magnitude increase in forward exchange reaction is required to see a separation of the integral responses.Initially, however, the two curves closely follow one another with the marginal diffusive advantage of Zn 2+ over Cu 2+ being insufficient to curtail the Cu 2+ binding dominance.At roughly twice the characteristic diffusion time, the two curves deviate markedly, with the CuPEI + integral curve increasing somewhat faster than linearly, while the ZnPEI + integral curve decreases linearly.The transition point is quite abrupt and common for the two curves.However, increasing k + 3 another three orders of magnitude results in complete separation of the two integral quantities, with the CuPEI + curve increasing monotonically with no apparent inection point and plateauing to a steady state value (of 0.622), well above the maxima reached by the other two cases (k + 3 values) during the simulation period.In contrast, the ZnPEI + curve reaches a maximum at approximately twice the characteristic time before decreasing linearly for another two characteristic time units to then asymptote to zero.The latter trend reects the plateauing of CuPEI + at the high positive value.We point out that the plateau reached by the CuPEI + complex content asymptotes to the theoretical limit (0.6306) of all binding sites being occupied by the one species; this is indicated by the upper black dashed line.The implication here is that the added effect of an additional exchange process drives the system to consider Cu 2+ being the preferred binding species.Unfortunately again, we are unable to track the time course of metal ion accumulation in the case of the smaller value of k + 3 ¼ 10 9 beyond the values shown.However, we anticipate that the two proles for this case will ultimately asymptote to the same steady state limits as for the k + 3 ¼ 10 12 proles, namely zero for ZnPEI + and 0.6306 for CuPEI + .
The detailed time and space development of the concentrations of the metal ion-ligand complexes, CuPEI + and ZnPEI + , for the case of k + 3 ¼ 1.0 Â 10 12 , clearly exhibits the more rapid increase in ZnPEI + concentration initially.However, as Cu 2+ spreads into the matrix and the concentration of CuPEI + increases, the depletion of ZnPEI + in favorable exchange with CuPEI + is apparent (data not shown).The exchange effect is complete in the long term (i.e., asymptotically) throughout the polymer matrix, but occurs sooner near to the edge of the Fig. 3 Integral graphs showing CuPEI + and ZnPEI + accumulation rates according to Scenario I. Line references as in Fig. 2.This data was generated assuming equivalent direct binding rate constants: and (practically) equivalent diffusion rates, while the different curves, in the order depicted by the arrows, result from increasing exchange reaction rate constants k + 3 ¼ 1.0 (pluses), 10 9 (circles) and 10 12 (diamonds), respectively.Other parameters are as given in Table 1.  1.
polymer layer as this region is rst to become saturated with Cu 2+ ions.
While we see a trend reecting the differing degree of binding, and specically the replacement of bound Zn 2+ for bound Cu 2+ , we have yet to encounter the observed event of initial dominant accumulation of Zn 2+ that is then overshadowed by the binding of Cu 2+ .If we now modify the diffusive process to disadvantage the diffusion of Cu by increasing the effective Cu 2+ radius to 4.5 Å, we nd, other things being equal, the qualitative behavior observed in the experiments of ref. 1.In Fig. 4, we compare the different cases of increased Cu 2+ radius R Cu ¼ 1.074, 4.5 and 6.5 Å, corresponding to a steady decrease in Cu 2+ diffusion constant.As in the case of Fig. 2, the rate of intrusion of Cu 2+ is slower as R Cu increases, allowing Zn to diffuse uninhibited into the matrix and bind with free PEI À sites.Only when a sufficient number of Cu 2+ ions have entered is there any decline in the rate of increase of bound Zn 2+ .However, in contrast to the results of Fig. 2, the rate at which the metal ions bind is dominated by Zn 2+ only initially, when diffusion of Cu 2+ is much slower (R Cu ¼ 4.5 Å compared with R Zn ¼ 1.069 Å).A cross-over point is reached beyond which the rate of uptake of Cu 2+ dominates over Zn 2+ binding through the action of the exchange mechanism (eqn ( 6)).Increasing R Cu to 6.5 Å reinforces the effect, although the peak value of ZnPEI + is only slightly higher and is reached at only a slightly later time (s ¼ 3.32 for R Cu ¼ 6.5 Å compared with s ¼ 3.19 for R Cu ¼ 4.5 Å).The cross-over point, however, is signicantly delayed in the slower diffusive case (s ¼ 5.26 for R Cu ¼ 6.5 Å compared with s ¼ 3.22 for R Cu ¼ 4.5 Å).Although we are not able to access later simulation times in the cases of R Cu ¼ 4.5 Å and R Cu ¼ 6.5 Å due to numerical difficulties, we again anticipate that the subsequent behavior of the proles will mimic the asymptotic behavior exhibited by the faster Cu 2+ diffusion case (R Cu ¼ 1.074 Å), since once again the only difference between the cases (difference in the diffusion constant of Cu 2+ ) only affects the rate of approach to steady state and not the steady state condition itself.
The combination of the observed short term behavior, depicted by the R Cu ¼ 4.5 Å and R Cu ¼ 6.5 Å results, with the expected asymptotic limit, demonstrated by the R Cu ¼ 1.074 Å result, correlates well with the behavior observed in the Lindén et al. experiments.

Scenario II
In the following calculations we shall assume the same total number density of absorption sites as assumed for Scenario I.This is to allow for some commonality between the two models.
In the majority of cases, however, half of these are assumed selective binding sites for Cu 2+ and half are selective binding sites for Zn 2+ .Again, in the majority of cases, the binding of metal ions will be limited only to their respective sites.As a point of model differentiation, this scenario models, more explicitly than in Scenario I, the situation of Cu 2+ binding having a direct impact on the binding of Zn 2+ through two mechanisms.
The rst mechanism is akin to the mechanism underpinning Scenario I, but rather than being considered an exchange process, the mechanism implies an induced conformational change to the polymer structure that destabilizes existing ZnPEI + complexes.A destabilized complex then undergoes a dissociation reaction to produce a free Zn 2+ ion and an inactive PEI À 2,ia binding site.The second mechanism is not considered at all in Scenario I. Here, the binding of Cu 2+ impacts directly on a bare PEI À 2,a binding site, converting it to its inactive state thus precluding any subsequent Zn 2+ binding.With regard to the former mechanism, we treat it here in two ways.First, by modelling explicitly the dissociation process as a separate, but coupled, reaction step, with its own dissociation rate constants (k AE equivalent to the rst, although it circumvents the need to identify possible k AE 4 values.In Fig. 5 we demonstrate the effect of increasing the effect of polymer reconguration triggered by the binding of Cu 2+ ions, manifested in the destabilisation of existing ZnPEI + 2,a complexes and the inactivation of bare PEI À 2,a sites.In this gure both treatments of the two-site model, with and without an intermediate dissociation step, are included to indicate their essential equivalence.Two points should be noted.First, with the separation into two site types, the theoretical limit of site saturation is equal for each binding species.This is indicated by the black dashed line and value of 0.3153.However, in contrast to Scenario I, as the binding of Zn 2+ decreases through the destabilisation mechanisms there is no corresponding increase in the level of Cu 2+ binding.Hence, in all the cases shown in this gure, the CuPEI + 1 integral curve is largely unaffected by the changes experienced by the ZnPEI + 2,a content in the matrix.The second point of note is that the curves were generated assuming equivalent ion diffusive rates (equal diffusion constants).Note also the absence of the experimentally observed behavior of initial dominance of bound Zn 2+ over bound Cu 2+ and subsequent reversal of binding trends.This sought-aer behavior, quite convincingly demonstrated in Fig. 6, is found when differences in diffusivities are introduced.Specically, as the rate of Cu 2+ intrusion into the matrix is successively reduced (by reducing the ion's diffusion constant relative to that of Zn 2+ ), there is a progressively stronger take up of Zn 2+ , at least in relative terms.However, the uptake of Zn 2+ never exceeds the theoretical limit (of 0.3153), and indeed it is always lower and slower than the amount taken up in the case of completely independent binding, which thus represents an upper limiting envelope in time for Zn 2+ uptake.This observation is to be contrasted with the results of Scenario I in which the ZnPEI + amounts can surpass the equal binding curve.Naturally, as the Cu 2+ diffusion constant is reduced, the ZnPEI + 2,a accumulation is larger, lying closer to the independent binding envelope.In all cases of Cu 2+ diffusion, however, the ZnPEI + 2,a curve attains a maximum aer which the level of bound Zn 2+ decreases linearly (at rst) before approaching the zero asymptote in the steady state limit, as more Cu 2+ binds to its specic binding site.The results in this gure were generated assuming an intermediate dissociation step.
It should be noted that as we have assumed an equal number of binding sites for Cu 2+ as for Zn 2+ , the steady state limit for the amount of ZnPEI + 2,a complexes will be zero.This limit would apply for any choice of underlying site densities satisfying the inequality n PEI 1,tot $ n PEI 2,tot .However, a nonzero steady state ZnPEI + 2,a limit exists under the current binding scenario should the underlying Cu 2+ binding site density n PEI 1,tot be less than n PEI 2,tot .Since experimental ndings indicate complete depletion of bound Zn 2+ , this would suggest the number of Cu 2+ binding sites to be greater than or equal to the number of Zn 2+ binding sites.
Fig. 7A shows the time and space development of the ZnPEI + 2,a complex under the condition of reduced Cu 2+ diffusion and signicant rates of inactivation and destabilisation.In Fig. 7B we show the corresponding increase in the density of inactive Zn 2+ binding sites, PEI À 2,ia .The latter quantity appears in response to the level of inux of Cu 2+ and thus the time and space dependent CuPEI + 1 density.For completeness, we show in Fig. 8 results of the model assuming no intermediate dissociation step in Scenario II, again for signicant inactivation and destabilisation reactions and for progressively reduced Cu 2+ diffusion.No major qualitative differences appear in comparison with the results shown in Fig. 6.

Summary remarks
With regard to Scenario I, we have shown that with the right combination of two inuential factors, a low relative Cu 2+ diffusion constant, a high Cu 2+ for Zn 2+ exchange rate, we can reproduce qualitatively, if not quantitatively, the experimental ndings of ref. 1. Naturally, it is necessary for the Cu 2+ binding rate/strength to be sufficiently high for the Cu 2+ binding to be at all competitive with Zn 2+ .However, the results indicate that it is not necessary for Cu 2+ binding to be signicantly stronger than the binding strength of Zn 2+ .The combination of the above two factors is sufficient to substantiate the hypothesis of a simple competitive mechanism (or mechanisms) for the reported experimental results.We remark that this scenario is the simpler of the two presented here.Advocates of Occam's razor would argue that this may be the most reasonable explanation.We remark also that nothing has  1.
actually been said of, nor does the model itself suggest, any specic physical or chemical processes that are being represented by the rate constants themselves.With regard to Scenario II, analogous summary comments can be made.A diffusion differential alone is not sufficient to explain the experimental ndings.Nor is a polymer matrix rearrangement triggered by Cu 2+ binding sufficient to account for the observed competitive binding behavior.The qualitative effect can be achieved through a suitable blend of these two mechanisms: signicant differential diffusion, and signicant destabilization of Zn 2+ -PEI À complexes and/or inactivation of Zn 2+ binding sites.
Incidentally, this is true whether Cu 2+ and Zn 2+ have different selective binding sites or they are attracted to and compete for the same binding sites (data not shown).
The features that distinguish the single site model from the two site model are points of differentiation that are potentially useful in aiding the interpretation of experiments on the real systems as well as for the design of new experiments for hypothesis testing.
In this work we have not made any attempt to accurately model the experimental system.Instead, we have taken the opportunity to explore the temporal and spatial details that physically underlie the competitive binding process.We have achieved the goal of identifying the essential processes that lead to preferential competitive binding.In our ongoing efforts, which will be the subject of future communications, we are undertaking the necessary rigorous quantitative comparison between theory and a more extensively characterized experimental system.In so doing, the model will be rened to improve its representation of the experimental system.This will include a more appropriate spatial account of ion diffusion through the porous matrix 29 and allowance for the possibility of matrix   1.
swelling during the absorption process.More accurate modelling, however, requires in turn a more accurate characterisation of the experimental situation as well as a more extensive survey of experimental conditions.Not only will this establish a solid benchmark of experimental data for a comparison, it will help to identify trademark behavior that could be used to distinguish between the different possible mechanisms responsible for the observed competitive ion absorption.
v 2 j vx 2 ðx; tÞ ¼ À rðx; tÞ 4p3 r ; 0\x\D; (26) in the polymer layer, and the Poisson-Boltzmann equation, N i expð À ez i jðx; tÞ=kTÞ; x .D; (27) external to the polymer layer.Here, n N i refers to the bulk ion density of mobile ion species i, k is Boltzmann's constant, 3 r is the relative dielectric permittivity of the matrix region and T is temperature in Kelvin.The two solutions are complemented by the boundary conditions of, rstly, zero surface charge at x ¼ 0 vj vx ðx; tÞ ¼ 0; at x ¼ 0; (28) secondly, the electroneutrality condition at innity, vj vx ðx; tÞ/0; as x/N;

Fig. 2
Fig. 2 Integral graphs showing the rates of accumulation of CuPEI + and ZnPEI + assuming independent binding reactions.Blue solid lines denote CuPEI + accumulation, while red broken lines depict ZnPEI + accumulation.Values shown are nondimensional.The thin black dashed lines indicate the level of saturated binding; the upper line denotes the level of saturated binding by one species, while the lower line indicates the level of equally shared saturated binding.The results were generated assuming direct binding rate constants of k+ 1 ¼ k + 2 ¼ 5 Â 10 11 and k À 1 ¼ k À 2 ¼ 4 Â 10 6and negligible exchange reaction, k AE 3 ¼ 1.0.In the order depicted by the arrows, we assume increasing effective Cu 2+ ion radius: R Cu ¼ 1.074 (pluses), 3.5 (circles) and 4.5 (diamonds) Å, respectively, or equivalently decreasing Cu diffusion constant: D Cu ¼ 2.26 Â 10 À6 , 6.94 Â 10 À7 and 5.40 Â 10 À7 cm 2 s À1 , respectively, compared with D Zn ¼ 2.27 Â 10 À6 cm 2 s À1 .Other parameters are as given in Table1.

Fig. 7
Fig. 7 Surface plots showing the temporal and spatial increase of bound species and inactivated Zn binding sites according to Scenario II with an intermediate dissociation step.The results assume the conditions of Fig. 6 for the extreme case of reduced Cu 2+ diffusion rate (D Cu ¼ 3.74 Â 10 À7 cm 2 s À1 compared with D Zn ¼ 2.27 Â 10 À6 cm 2 s À1 ).The subplots show (A) ZnPEI + 2,a density and (B) PEI À 2,ia density.All variables are nondimensional.
and Zn 2+ under concentrated electrolyte conditions, are not too dissimilar to what can be deduced from the ideal Stokes-Einstein relation assuming h is the viscosity of water and where R is the van der Waals radius of the ion.Taking van der Waals radii of 1.074 Å, 1.069 Å and 1.18 Å for Cu 2+ , Zn 2+ and Cl À gives the values of 2.26 Â 10 À6 , 2.27 Â 10 À6 and 2.06 Â 10 À6 cm 2 s À1 for the ions' respective diffusion constants.As we have no knowledge yet of the effective diffusion constants of these ions in the polymer matrix, we shall utilize eqn (40) to derive diffusion constants based on effective ion radii, in order to readily draw a comparison with diffusion in the bulk solution.