Lability and mobility effects on mixtures of ligands under steady-state conditions

Analytical solutions for the steady-state flux arriving at an active surface from a mixture (in which one active species reacts with non-active ligands in the medium) can be helpful in a variety of problems: voltammetric techniques, heterogeneous processes in reactors, toxic or nutrient uptake, techniques of diffusive gradients in thin films (DGT), etc. Under any geometry that sustains steady-state, a convenient combination of the reaction-diffusion equations leads to a simpler formulation of the problem for arbitrary diffusivities of the species and arbitrary rate constants of the first order conversion between the active species and the non-active species. The resulting problem can be characterised in terms of a list of dimensionless parameters involving the kinetic and mobility properties of each species. A lability degree for each 1:1 complex in terms of the surface concentrations leads to: i) a lability criterion specific for each complex in the mixture and ii) the assessment of the relative contribution of each complex to the resulting flux. Semi-infinite spherical diffusion (as in the Gel Integrated MicroElectrode, GIME, biouptake modelling of micro-organisms, etc.) is specifically considered and some consequences of its full analytical solutions are discussed. Published in Physical Chemistry Chemical Physics 2003, vol 5, p 5091-5100 DOI: 10.1039/b306172h reprints to galceran@quimica.udl.cat


Introduction
Diffusion coupled with chemical reaction in the solution appears in a variety of problems in different scientific areas such as: i) Cellular nutrition [1][2][3][4] and biouptake of pollutants by micro-organisms [5][6][7] .In many cases, the molecules of nutrients or toxins can be complexed by ligands present in the medium while the diffusive transport occurs.The bioavailability of the substance of interest will depend not only on the nature of the active surface (the cellular membrane in this case) and the dynamics of diffusion, but also on the extension and the kinetics of the complexation processes in the medium (i.e. its kinetic speciation [8;9] ).
ii) Heterogeneous catalysis and electrode reactions.The problem is also relevant to Chemical Engineering [10] because the efficiency of the heterogeneous reaction is affected by diffusion and reaction in the solution phase.
In all of these systems, the responses due to the arrival of the active substance at the surface depend not only on the equilibrium values of its complexation with the ligands (thermodynamic speciation), but also on the particular values of the various association and dissociation rate constants, i.e. on the dynamic speciation.Two limiting cases are easily recognised: the inert limit (no dissociation of the complexes) and the fully labile limit (the equilibrium relationship holds at any point of the solution).Thus, the term lability can be understood as a continuous property (between the aforementioned limits, passing through different amounts of partial lability) measuring the contribution of the dissociation of the complexes to the response flux [25;26] .In this context, lability depends not only on solution properties (values of diffusion coefficients, complexation and rate constants, etc.) but also on the time scale of the experiment, the geometry of the surface, the thickness of the diffusion domain, etc.
The aim of this work is to study steady-state fluxes arriving at an active surface due to diffusion and reaction of an active species and non-active ligands, taking into account all the dynamic properties of the mixture, including different mobilities and no restriction for the values of the association/dissociation rate constants.We extend the procedure of Turner and Whitfield [12] who developed a very general solution of multicomponent reaction, convection and diffusion (with all components sharing a common diffusion coefficient).By considering diffusion as the only transport phenomenon, we can deal with unequal mobilities (i.e.unequal diffusion coefficients) for the different complexes, so that macromolecular complexation -which is essential for understanding environmental problems [27] -can be realistically modelled.
For simplicity, but without losing the generality of the above mentioned applications, we use the voltammetric terminology.The active surface is the electrode, which is the sink of the metal ion M (active species).In the solution, M can be complexed by h ligands to form h different complexes, none of which is (electro)active on the electrode surface.In the given examples, we will restrict ourselves to fluxes in diffusion limited conditions (null concentration of M at the active surface).The extension of the treatment here expound to the case of a fixed arbitrary value for the concentration of the active species at the active surface (due to phenomena such as Michaelis-Menten internalisation or a Nernstian relationship between species at the surface) is straightforward and is not detailed here for simplicity reasons.
In the first section of this work, we lay down the re-formulated dimensionless equations that should be solved under any specified geometry.In section 2, we introduce the concept of degree of lability for each complex, which allows the formulation of lability criterion for each complex.In section 3, we concentrate on spherical geometry with semi-infinite diffusion, presenting closed expressions for a mixture of two ligands.
Derivation of expressions for one dimensional geometries are included in Appendix A, while the symbols used are gathered in Appendix B.
The results of this article can be applied to discrete mixtures of ligands, either arising from the mixing of a number of different real species or from the consideration of sites in macromolecules as formal species which can have different affinities (e.g. if different functional groups are involved in different kinds of sites) but have the same mobility.

1.1.-Setting of the problem
Let a multicomponent mixture consist of one active species (the free metal ion M in the adopted terminology) and h different ligands designated L 1 , L 2 ..L h , exhibiting h simultaneous complexation reactions of the type a, d,

M+L ML
where subscripts a and d indicate association and dissociation respectively.See Fig 1 for an outline of the kind of system considered.Assuming ideal dilute behaviour, the equilibrium relationship (fulfilled in the "bulk", denoted with superscript *, either at a finite or infinite distance from the active surface) is: Assume that the ligand concentrations are large in comparison with that of M so that their value can be approximated as constants * L i c at any point of the solution (i.e.ligand excess conditions).The excess of ligand condition for a given ligand i is clearly fulfilled ; but the ligand of excess approximation can yield good results for the flux even at lower ligand-to-metal ratios.When ligand excess conditions apply, the kinetics of interconversion of M and ML i follow a pseudo first-order law and the problem becomes linear.
Assume a steady-state mass-conservation equation for each complex ML i with diffusion as the only relevant transport phenomenon: where the (dimensional) Laplacian operator depends on the particular geometry considered.For the free metal ion, the mass-conservation equation reads: In order to work preferentially with dimensionless parameters we define: which is an index of the mobility of the particular complex i with respect to the metal ion (active species).For the association kinetics, we introduce where ℓ stands for some characteristic length of the active area (such as the radius of the sphere or of the disc electrode ) [28;29] or a characteristic length of the system (such as a diffusion layer thickness).κ a,i is a dimensionless association parameter [18;30] , Damköhler number or Thiele modulus [10;31-34] , comparing the diffusion efficiency of the metal and its association rate constant with a given ligand L i .One noticeable fact is that κ a,i depends on the characteristic length: for decreasing values of ℓ , κ a,i decreases indicating that the role of diffusion increases with respect to the role of kinetics [18] .For a fixed metal and a fixed ℓ , the κ a,i -values are just proportional to the corresponding k a,i of the reaction of M with the different ligands.
For notation simplicity we also define 0 a, =1 Thus κ 0 indicates the added effect of all the ligands.
Similarly, for the dissociation process which can be seen as the dimensionless dissociation parameter comparing, for a given complex, the effectiveness of the processes of its dissociation and its diffusion.
Notice that a, ' d, We also introduce the normalised concentrations: where now the Laplacian operator is dimensionless.
We consider 2 response functions: i) the gradient of normalised concentration of metal grad(c 0 ) at the electrode surface and ii) the normalised flux of metal φ, which is the ratio between the total flux measured with and without ligand [35;36] for the same total metal which in voltammetry can be simply obtained as the ratio of currents with and without the ligand presence [15] .

1.2.-Method of solution for any geometry
The set of h+1 differential equations ( 13) and ( 12) can be written in matrix form as: We seek a combination of equations ( 12) and ( 13) so that they are re-written as where f j are some combinations of the concentrations c i ; i.e. the goal is to uncouple the differential equations for some formal species concentrations f j (also called "quasispecies" in the theory of macromolecular evolution [37] ).In matrix form: (18)   where D is a diagonal matrix.If A can be rendered diagonal, there exists a matrix S so that: A S D S (19)  Thus, the sought relationship between the concentrations can be expressed: or Finding the combination The elements of the diagonal matrix D are the eigenvalues of the matrix A, which can be found as the zeros of the characteristic polynomial of A. Following the usual procedure [38] , we have to solve: (22)   where I is the identity matrix.In the determinant, we add all the rows to the first one taking into account the definition (7); we factor out t from the first row and factor out ( ) − from each row j (2≤j≤h), so that eqn.( 22) can be expressed as ( ) ( ) which, after subtracting all the rows from the first one, leads to One solution of the previous eqn. is t=0, which can be labelled as the eigenvalue n 0 .The other h eigenvalues n 1 , n 2 ,..., n h required in eqn.( 17) can be found from a, which is a polynomial equation of degree h in t.We stress that n 1 , n 2 ,..., n h are function of κ d,i and κ d,i , but totally independent of the geometry.
It is worth noticing that all the values n 1 , n 2 ,..., n h are positive, each one lying between two κ d,i : This property can be deduced from considering the graphical solution of eqn.(25) in the following form a, 1 d, each n j is the intersection between the bisectant of the first quadrant (r.h.s. in eqn.( 27)) and a sum of hyperbola-like functions (l.h.s. in eqn.( 27)), as seen in Fig 2 for the simple case of h=3.The knowledge that there is one eigenvalue n j between two consecutive κ d,i can be exploited for accurately determining the n j -values.On the other hand, since all the eigenvalues are different, matrix A can be diagonalized and all the construction is proved consistent.
In general, each n j depends on the whole set of κ a,i -values and κ d,i -values, as indicated in eqn.(25).For a given formal species with concentration f j , the parameter n j can be physically understood, according to eqn.(17), as the non-dimensional kinetic constant of disappearance of that formal species.Larger n j -values correspond to formal species having thinner reaction layers.
Suppose that we know (e.g.numerically) the eigenvalues.Now, we turn to the problem of finding the coefficients s ij in S. As the columns in matrix S are the eigenvectors [38] , we just need to solve the eigenvalue problem: ( ) (28)   which is the system As the previous system is undetermined, the first equation can be skipped and, as the eigenvectors must be non-null, we can choose 0 1 0 which leads to 0 0 and the other coefficients can be easily found from the remaining equations in (29):

Computing the flux
We assume that M is the only active species on the active surface (electrode).Thus, A useful result can be obtained by adding up where eqn.( 25) has been used.Thus, ( ) ( ) which indicates that the flux of M can be easily computed once f 0 is known.

Outline of the solution procedure
Once the roots n 1 ,..., n h in ( 25) are found, one must solve (17) with the appropriate boundary conditions: BOUNDARY CONDITIONS FOR f j AT THE ACTIVE SURFACE: On one hand, at the electrode surface, we have no flux given by ( 33) corresponding to c i (i>0).These conditions, referred to f i via (21), become: ( ) 0 grad( ) 0 grad 0 , active surface If we are under diffusion limited conditions, then c 0 =0, which referred to f j via (31) becomes 0 0 0 active surface else some value should be provided for c 0 (directly or from complementary information such as Nernstian equilibrium).Notice that the boundary conditions to apply to each f i at the active surface are implicit within the summations ( 36) and (37), i.e. the combination (changing from c j to f j ) that simplifies the differential equations mixes up the boundary conditions everywhere (and thus at the active surface, here considered).So, in order to prescribe the boundary conditions on a given f j one needs, in general, some kind of inversion of the matrix S (see eqn. (51) in ref. [12] ).
BOUNDARY CONDITIONS FOR f j IN THE BULK: On the other hand, more boundary conditions are needed, even for 1-D geometries.In some cases, fixed bulk conditions are prescribed, either at a finite or infinite distance.In these cases, eqn.(34)   leads to the boundary condition where definitions (10) and ( 11) have been applied.
For the other formal species, we notice that eqn.( 21) under bulk conditions, reads once ( 11), ( 30) and ( 38) have been taken into account.The previous system of equations (39) is homogeneous leading to the trivial solution * 0 1 which provides the boundary conditions to prescribe.

The parameters defining the problem
The response given through grad(c 0 ) (but not through φ ) corresponding to equations ( 3)-( 4) with boundary conditions ( 36)-( 40) could be defined with 3 lists of parameters: κ a,i , κ d,i and ' i i K ε , apart from geometrical factors linked to the boundary conditions such as those related to the diffusion domain not included in ℓ .Recalling the relationship between the association and dissociation parameters given by ( 9), we notice that only 2 (of the original 3) lists are necessary.

2.-One-dimensional geometries: the degree of lability
As mentioned in the introduction, the concept of lability is related to the availability of the active species arising from the dissociating complexes.In order to qualitatively assess such contribution, the concept of degree of lability for the particular case of only one complex under diffusion limited conditions in spherical geometry was introduced [18] .The degree of lability for that case was shown to be directly related to the normalised concentration of complex at the active surface: where superscript 0 indicates the concentration at the active surface.Indeed, the normalised concentration of the complex species at the active surface is a good measure of the capacity of the complex to dissociate when c M is depleted and thus of the lability of this complex (for planar geometry results, see Figs 1, 2 and 3 in ref. [19] ).
The gradient of metal concentration at the origin for this spherical case with only one complex could have been written as: ( ) We now proceed to generalise this concept to a mixture of ligands in one-dimensional geometries.The first step consists in recognising the key role of the fraction of each complex normalised concentration at the active surface on the flux (and, so, on φ).
Adding up the h eqns. in (12) with eqn.(13), one obtains, As shown in appendix A, restriction to any one dimensional geometry that sustains steady state (such as planar, cylindrical or spherical diffusion in a finite domain or spherical diffusion in a semi-infinite domain) allows a simple solution of eqn.(43) and the flux to be written in terms of the differences of complex concentrations between the bulk and the active surface From the definitions, these differences can be re-organised as: ( ) It can be shown (see eqns. By comparing eqn.( 45) with (42), we suggest the extension of the concept of degree of lability [18;25] to each complex of the mixture under any one-dimensional geometry and any steady state regime (i.e.not necessarily diffusion limited conditions) as: However, for simplicity, from now on we only consider here the diffusion limited case In order to illustrate the suggested ideas, we consider a model mixture whose parameters are gathered in Table 1.There are two complexes: one (ML 1 , a macromolecular species) with strong stability constant, low abundance and low mobility and the other (ML 2 , formed with a small ligand) with weak stability constant, large abundance and same mobility as the metal ion.For these conditions, with a spherical active surface of radius 10 µm, the degree of lability of the complexes (computed with expressions given below in section 3) is ξ 1 =0.06 and ξ 2 =0.33 (see Table 2), as indicated by the height of the corresponding bars in Fig 3.
Using (45) and (46) in eqn.( 14), the normalised flux (under diffusion limited conditions) for any one-dimensional geometry can be written To apply this expression numerically one still needs to solve the problem as outlined in section 1, but expression (47) (or ( 45)) indicates that the potential contribution of each complex ( ' i i K ε ) to the overall flux is modulated by the value (relative to the bulk) of its concentration at the active surface (see definition (46)).Indeed, if we plot a bar for each complex i at the abscissa  2 for the numerical values), one concludes that the contribution to the flux of complex ML 1 is much larger than the contribution of ML 2 , despite ML 1 has a lower degree of lability.In this example, the contribution of the free metal to the overall flux (equivalent to the area between the bar at ' 1 i i K ε = and the ordinate axis) is negligible in comparison with the contribution of the complex.
A further conclusion from eqn. ( 47) is that the potential contribution of a particular complex depends on ε << , the presence (and lability) of the complex i has a negligible impact on the total flux because the metal has a larger weight, so the assumption ' 1 i i K ε >> can be used whenever a simplification is necessary.
Another application of the degree of lability is the possibility of extending the lability criterion concept [14;18;25;39;40] to each complex in the mixture.One can impose as the condition a specific complex ML i has to fulfil in order to be predominantly labile.

3.1.-Solution of the problem
In this case a relevant characteristic distance is the electrode radius ( 0 r = ℓ ).Thus, a normalised co-ordinate can be defined as: The solution of eqn.(17) with j=0 (i.e.n j =0), imposing only the bulk boundary condition (38) and leaving an unknown coefficient (T 0 ) to be determined from the other boundary condition, is: while the non-flux condition for the complexes -see (36)-becomes Combining eqn.( 52) and ( 53), one obtains d, which is a linear system of equations allowing T 1 ,..., T h to be found.The sought gradient can be computed as After some algebra combining eqns.( 10), ( 21), (32), n 0 =0, (50), (51), ( 52) and ( 9), one finds that the degree of lability for the spherical active surface with semi-infinite diffusion can be computed as:

3.2.-Studying the case of two ligands
With the previous formalism, the expression for the flux found for h=1 can be retrieved [18] .We consider, then, the next case corresponding to h=2.
From (25), the general expressions for n 1 and n 2 (regardless of the geometry) are where we have chosen the sign minus in front of the root for n 1 and the sign plus for n 2 , so that n 1 < n 2 (to be consistent with ( 26)).
Solving (54) with h=2 and replacing into (56), we obtain the general expression for the dimensionless flux: As a first particular case, consider that the complex ML 2 is totally inert (κ a,2 , κ d,2 and n 1 tend to zero while n 2 tends to κ a,1 +κ d,1 ) while there is no restriction on the lability of ML 1 .Then, ( ) A second particular case consists in the complex ML 2 being fully labile (κ a,2 and κ d,2 tend to infinity) while there is no restriction on the lability of ML 1 .Then, It has been demonstrated [18] that decreasing r 0 results in the loss of lability in systems with only one complex, due to diffusion enhancement for smaller r 0 .The same conclusion can be reached for mixtures of ligands, as of both complexes at this radius can also be seen in the normalised profiles depicted in Fig 5, where we also notice that the reaction layer is much thinner than the diffusion layer [19] : the former is seen to disappear around r/r 0 ≈ 1.00002 in the figure, while the latter extends beyond the plotted region.
For some combination of the parameters (see Fig 6), an almost successive loss of lability of each complex can be seen in the chosen logarithmic representation.This fact indicates a potential use of electrodes with variable radii: the kinetic characteristics of the complexes could be easily retrieved from the succession of recorded steps in the measured fluxes.
For prepared mixtures (i.e.not arising from ligands having a fixed proportion of different types of sites in the same molecule), one has experimental control on individual * L i c , which results in a capability of changing ' i K and κ a,i which determine a variation in the measured fluxes.Let us assume that we prepare several mixtures with 2 kinds of ligands having each ligand a different complexation site, so that the summation of the concentrations of both sites is fixed in all the mixtures (i.e. ).We caution that experiments done in this way must fulfil the excess of ligand condition for both ligands at each prepared mixture.We see the appearance of a maximum of the flux for a certain intermediate molar fraction ( ) where the kinetic constants of both ligands are not very different, the more mobile ligand is less labile ( ε 1 >ε 2 ; k d,1 <k d,2 ) and we have chosen a suitable radius.From left to right in Fig 7 we increase the proportion of L 1 which diffuses faster than L 2 , but dissociates slowlier than L 2 .We could interpret the figure saying that adding L 1 increases φ at the first additions because of its larger ε 1 , but decreases φ at further additions because of its smaller k d,1 .
We can combine the 2 experimentally controlled variables: r 0 and composition, as seen in Fig 8 .The continuous line here is the same as in the previous figure 7, but now looks flatter because of the larger axis scale.If we increase r 0 (dashed line), the complexes are more labile and the limiting step is diffusion: so increasing L 1 yields larger φ-values for most of the composition range.If we decrease r 0 (dotted line), diffusion is no longer the limiting step, so the addition of L 1 only means a less labile mixture, and, thus, φ decreases.
Appendix A: The flux depends on the differences between bulk and active surface concentrations.
For any one dimensional geometry, equation ( 43) can be written [41] 2 where γ is a convergence parameter taking the value 1, ½ or 0 for spherical, cylindrical or planar geometry respectively and ρ is the distance from the centre of symmetry of active surface normalised with respect to the characteristic distance ℓ .The solution of eqn.(A-1) is [41] : where the constants B 1 and B 2 are to be determined from the boundary conditions.We prescribe fixed concentrations at the active surface ρ 0 : With these boundaries conditions -and recalling the absence of gradient for the complexes as expressed by (33)-one computes the gradient of c 0 at the active surface as: ( ) which mean that the gradients of M are proportional to the differences between bulk and surface concentrations of all species.
Additionally, as eqn.( 34) can be re-written as the solution of the Laplace eqn.(17) for i=0 (n 0 =0), i.e. f 0 , can be obtained by simply dividing eqn (A-2) (or (A-3)) by ; 0 0 11)   where I and K are first and second kind modified Bessel functions.
Then, the gradients can be found to be:    (27) showing that all the eigenvalues are non-negative.
The dotted line stands for the r.h.s of eqn.(27), while the continuous line stands for its l.h.s.Dashed lines stand for the asymptotes.shows the contribution of this complex to the total flux according to eqn.(47).
Parameters from Table 1 and r 0 =10 µm.For completeness purposes, the contribution of the metal (in this case negligible area where there is no room for filling) is also included.

Metal
given by the degree of lability for this complex, then the area of the rectangle determined by this bar and the axes (see different fillings of areas in Fig3) yields the contribution of the chosen complex to the flux indicated by eqn.(47).In order to reproduce the leading 1 in the numerator of eqn (47), a virtual contribution of the free metal can also be included in the figure by adding a bar of unit height at ' solution of eqn.(17) taking into account the bulk boundary condition (40) can be written (e.g.see Appendix A): impose the boundary conditions on the electrode surface (ρ=1).In case of diffusion limited flux of M -see (37)- Fig 4 shows: the normalised current drops when the radius shrinks due to the diminishing of the complex dissociation contributions to the measured flux.In Fig 4 we see a dramatic change in the dynamic behaviour of this particular system for radii around 1 mm.The partial lability and at a certain distance (ρ=ρ*, finite or infinite) to formulate the solution in terms of the differences of those concentrations (even though the values of 0 i c are a priori unknown).

Fig 1 :
Fig 1 : Outline of the multicomponent diffusion and reaction (see scheme (1)) towards a

Fig 2 :
Fig 2: Graphical solution of eqn.(27) showing that all the eigenvalues are non-negative.

Fig 3 :
Fig 3 : Plot of the degree of lability ξ i vs. the product ' i i K ε

Fig 4 :
Fig 4 : Impact of changing the electrode radius on the lability.Continuous line stands

Fig 6 :
Fig 6 : The predominant loss of lability of each of the 2 complexes of the mixture as r 0 decreases (continuous line) leads to an almost two-step wave.Parameters: k a,1 = 10 6 mol -

Fig 7 : 1 m 3 Fig 8 :
Fig 7 : Variation of the normalised flux φ with the molar fraction -x 1 -of ligand L 1 in a binary mixture leading to the appearance of a maximum.Parameters k a,1 =k a,2 = 10 2 mol -1