The Becker-Döring equations with monomer input , competition and inhibition

We investigate the Becker-Döring model of nucleation with three generalisations; an input of monomer, an input of inhibitor and finally, we allow the monomers to form two morphologies of cluster. We assume size-independent aggregation and fragmentation rates. Initially we consider the problem of constant monomer input and determine the steady-state solution approached in the large-time limit, and the manner in which it is approached. Secondly, in addition to a constant input of monomer we allow a constant input of inhibitor, which prevents clusters growing any larger and this removes them from the kinetics of the process; the inhibitor is consumed in the action of poisoning a cluster. We determine a critical ratio of poison to monomer input below which the cluster concentrations tend to a non-zero steady-state solution and the poison concentration tends to a finite value. Above the critical input ratio, the concentrations of all cluster sizes tend to zero and the poison concentration grows without limit. In both cases the solution in the large-time limit is determined. Finally we consider a model where monomers form two morphologies, but the inhibitor only acts on one morphology. Four cases are identified, depending on the relative poison to monomer input rates and the relative thermodynamic stability. In each case we determine the final cluster distribution and poison concentration. We find that poisoning the less stable cluster type can have a significant impact on the structure of the more stable cluster distribution; a counter-intuitive result. All results are shown to agree with numerical simulation.


Introduction
In 1935 Becker & Döring [1] presented an enduring model of nucleation where clusters form by the addition, or subtraction, of single particles (monomers) with no interaction between larger clusters.Such larger clusters evolve by maintaining a dynamic balance of monomer aggregation and fragmentation.Modelling this process as a series of chemical reactions and denoting an r-sized cluster by X r , we have For each reaction there are two reaction rates to prescribe, we denote the forward rate by a r and the reverse by b r+1 , both non-negative.Defining J r as the net flux from cluster size r to r + 1 and x r (t) as the concentration of clusters X r at time t, we express the system by where originally x 0 (t) was assumed to be such that ẋ1 = 0 [1].Later Penrose [2] generalised this by setting x 0 (t) = 0, thus ensuring the conservation of density These modified equations are still referred to as the Becker-Döring equations.For certain aggregation and fragmentation rates the existence and uniqueness of a solution to (2)-( 4) have been demonstrated by Ball, Carr & Penrose [3] for densities below a critical value; furthermore this result was subsequently generalised to arbitrary initial data by Ball & Carr [4].The asymptotic solution for a variety of aggregation and fragmentation rates has been described by Wattis & King [5], King & Wattis [6] and Wattis et al. [7].
Various aspects of the Becker-Döring equations have been investigated, including the existence of metastable solutions by Penrose [8], the aggregation-dominated regime by Carr [9] and the difficulties in numerically modelling metastable systems by Carr, Duncan & Walshaw [10] and Duncan & Soheili [11].The self-replication of micelles and vesicles, including the size-templating matrix effect, have been successfully modelled by novel generalisations of the Becker-Döring model [12,13,14,15].Additionally the origin of RNA has been studied by Coveney & Wattis [16].While being widely applicable, the Becker-Döring equations make the restrictive assumption that only monomers may interact with clusters.Smoluchowski [17] proposed a more general model allowing all cluster sizes to aggregate, and for a cluster to split into uneven fragments.Blackman & Marshall [18] exploit the Smoluchowski equations to study scaling behaviour in essentially the Becker-Döring regime.Da Costa [19] generalised the Becker-Döring model to include dimer interactions; a model subsequently analysed by Bolton & Wattis [20].Thus the Becker-Döring model with either x 0 (t) = 0 or x 0 (t) such that ẋ1 = 0 has been extensively studied and generalised.We propose to investigate a more general system which maintains a constant influx of monomers, that is x 0 (t) = x 0 , which is relevant to many industrial processes which rely upon continuous flow reactions rather than production in batches.
A constant input of mass can be balanced either by removing mass at larger cluster sizes or by adding a poison influx; and we model the latter.The inhibitor is consumed in the action of poisoning a growing cluster; once poisoned, a cluster has no further interactions and is assumed to be removed from the system.Previously inhibition has been studied by Wattis & Coveney [21], where clusters above a chosen size could be neutralised by an inhibitor.However, this description lacked an influx of poison, and monomer, and rather considered the case where the system initially has a finite stock of monomer and poison.
Finally we generalise the above poisoning model by allowing competition for the monomers which can form clusters of two morphologies and allow only one morphology to be poisoned.Preliminary work on this system has been reported by Wattis [22] where the constant density system (x 0 (t) = 0) was analysed and included the possibility of two morphology of clusters developing from the monomers, the first is denoted by X r and the second Y r .Thus the reactions that were allowed had the form No cluster can change morphology from X to Y or vice-versa; the only way mass can change from one form to the other is by the stepwise break-up of one cluster entirely into monomers (which have no morphology) and the subsequent re-aggregation of monomers in the other form.The system analysed in [22] was shown to have the same basic properties as the original Becker-Döring system, that is a conserved density, a unique equilibrium solution and a Lyapunov function.Again this model lacks the influx of monomers and poison which we will consider here.This introduction concludes with a description of the general model of monomer input, competition and inhibition, and how this model is truncated so that numerical simulations may be performed.Due to the complexity of the model we study each generalisation to the original Becker-Döring system in turn.In Section 2 we formulate the model which includes a constant influx of monomers, and assumes that there is no poison and only one morphology of cluster can form.We derive the solution that is approached in the large-time limit and calculate the large-time kinetics.In Section 3 we add a poison influx to the model, still allowing only one morphology of cluster.We identify a critical level of poison influx and solve the system when the poison addition rate is above and below this threshold.Finally, in Section 4, we consider the full model, with two morphologies forming, a poison influx which affects only one type of cluster and a monomer influx.In particular we consider how poisoning the less stable cluster type can influence the concentration profile of the more stable cluster.The paper concludes with a discussion of the results in Section 5.

General model
In the format of chemical reactions if we permit the X r clusters to be poisoned, we include as well as the reactions (6).The poisoned clusters, U , are assumed to be completely inert and have no further interaction with the system and so are ignored in the derivation of the kinetic equations.We define the rate at which clusters are poisoned to be k r and applying the law of mass action to ( 6) and ( 7) we obtain the infinite set of differential where p 0 (t) is the input rate of the poison and x 0 (t) is the input rate of the monomer.The concentrations of X-clusters are denoted x r (t) and those of the Y morphology by y r (t).This system is much more complex than the original Becker-Döring model, with the addition of not only competition but also an influx of monomers and poison.To simplify the system we assume throughout this paper that all the aggregation, fragmentation and poisoning rates are size-independent, that is a r = a, b r = b, α r = α, β r = β and k r = k; also that the influx of monomers and poison are time-independent, so that x 0 (t) = x 0 and p 0 (t) = p 0 .In addition we assume that the initial conditions are x r = y r = p = 0 ∀r, that is there are no clusters, monomer or poison present initially.
We study several simplifications of ( 8)-( 13) before finally including all the terms in Section 4.

Numerical Simulation
To solve the system numerically we truncate the system at a finite size r = N .However, due to the constant monomer influx, we expect a steady-state, or a borderline equilibrium, solution to arise (x r = x 1 or y r = x 1 ∀r).We assume that J N = ax 1 x N −bx N and I N = αx 1 y N − βy N ; this ensures that if x N , y N = 0 then J N , I N = 0 and hence the equilibrium solution will be correctly reproduced and also if x r = x 1 or y r = x 1 then this steady-state, or borderline equilibrium, solution will also be correctly reproduced.With these boundary conditions we numerically solve We solve the truncated system of equations by use of ode23s and ode15s of the Matlab 6 package [23]; these are special solvers for stiff systems which are accurate to large times, as required by some of our simulations.All numerical simulations referred to in this paper are of the above form with the relevant aggregation, fragmentation and poisoning rates and influxes as given by the particular example.

Constant flux
In this Section we study the Becker-Döring model with a constant influx of monomers; however, in this simplification we do not allow for a poison influx or for competition for the monomers and so the equations to be studied are Assuming that no clusters are present initially, the dynamics arise due to the steady influx of monomers.In the absence of input the equilibrium solution is given by where θ = ax 1 /b which has density ρ = x 1 /(1 − θ).Since the density in the system (20)-( 21) satisfies ρ = x 0 t we expect that θ −→ 1 as t −→ ∞.Numerically solving the system reveals that as more mass is injected into the system x 1 rises until θ approaches unity and the borderline equilibrium solution x eq r = x 1 = b/a is approached, as mass aggregates to increasingly large cluster sizes.
Consider the zeroth moment, defined by which, with the equations for ẋr (21) for r = 1, 2, . .., yields Obviously this is dependent on the monomer concentration but we know that the monomer concentration is tending to x eq 1 = b/a and so we assume that x 1 (t) = x eq 1 +q 1 (t) with q 1 (t) −→ 0 as t −→ ∞, and hence obtain θ = 1 + aq 1 /b which in equation ( 24) yields, at leading order, and leads to the scalings To consider the correction terms of equation ( 24) we need to include more terms in the expansion of M 0 and x 1 , so that and so we obtain To balance all the terms in this equation, we require that µ = 1/2 and ν = 1/2; numerical simulations confirm that µ = 1/2.We proceed by investigating the evolution of the size-distribution x r (t) as t −→ ∞.We assume that x 1 = b/a + q 1 t −1/2 and x r = bψ(r, t)/a with ψ(r, t) −→ 1 as t −→ ∞, where ψ = 1 for r = 1 and ψ −→ 0 as r −→ ∞.From equation (21), for ẋr , we obtain ( This equation has a self-similar solution which can be written in terms of the selfsimilar variable η = r/ √ t, yielding the solution ψ = A erfc r/2 √ bt − s 1 , where b and we determine A from the boundary condition ψ = 1 at r = 1; that is A = 1/erfc(−s 1 ), and so finally In summary, we have considered a simplified system which neglects competition and poisoning effects, but includes an influx of monomers.Numerically we found that at small times the system tends to a local equilibrium solution, and over large times this solution tends to x r = x 1 , with increasingly large clusters being produced.The timescale over which this occurs has been investigated and it was found that the monomer concentration tended to b/a with correction term decaying with O(t −1/2 ), and that x r is given by the similarity solution (30) as t −→ ∞.

Inhibition
Having investigated the model without an inhibitor we proceed to consider the more general case of a non-zero influx of poison, p 0 > 0, but maintain only a single morphology of cluster (i.e.no competition for monomers).Hence we consider the equations We seek a solution to this system of equations which balances the monomer influx by the poison influx.The behaviour depends on the relative value of the monomer to poison influx and all possible ratios are considered below.

Case
First we consider the parameter regime 2p 0 < x 0 ; that is, a relatively small influx of inhibitor.The addition of poison will permit the existence of steady-state solutions in the large time limit and for these solutions we assume that x eq r = θr−1 x x 1 and the effect of poisoning suggests that θx < 1. Requiring ẋr = 0 with (32) yields which is solved by Solving ṗ = 0 in (33) gives and finally if ẋ1 = 0 then (31) implies that and these three equations are solved to give θ, θx and p. Substituting (36) into (34) implies and we discount the possibility that in general θx = 1 since this leads to a contradiction in (34); thus θ − θx = ap 0 /b 2 θ and we obtain Given that x 0 > 2p 0 , equation (35) yields with x 1 = bθ/a; we take the positive root of (40) so that θ −→ 1 in the limit k −→ 0. Thus the final value of the poison is We test this analytical result by comparison with a numerical simulation.In the large-time limit formula (41) is correct as demonstrated in Figure 1A where p(t) − p final is plotted against time, where p(t) is the numerical result; this difference tends to zero as t −→ ∞ and we note that, with a log scale on the vertical axis, the graph is linear, indicating exponential decay.In Figure 1B we compare the numerical and analytical solutions and find they agree with good accuracy; furthermore we find numerically that the concentrations x r decay exponentially to the given steady-state solution x r = x 1 θr−1 x .Given the condition x 0 > 2p 0 then this analysis holds, and we note that θx is independent of a and b (numerically confirmed); however, x 1 is dependent on a and b.If p 0 = 0 then θx = θ = 1 as t −→ ∞ as previously obtained.We also note that the average cluster size, M 1 /M 0 = 1/(1 − θx ) is independent of the susceptibility of poisoning (k), but the final poison concentration is proportional to 1/k.From (40) we see that θ −→ 0 as p 0 −→ x 0 /2 and (41) implies p final −→ ∞ as p 0 −→ x 0 /2; we will return to these results when we consider the case x 0 = 2p 0 .

Case B:
Increasing the level of poison will have a non-trivial impact on the system.Numerical simulations suggest that all concentrations tend to zero, except for x 1 (t) which tends to a constant, and p(t) which grows without limit in a linear fashion.To investigate the large-time kinetics we make the assumption that x r (t) ∼ z r /t γr and p(t) = p 1 t; further assuming that γ 1 = 0. Thus equation (32) becomes as t −→ ∞ the dominant terms are the first and last term on the right hand side and by balancing them we find that γ r − 1 = γ r−1 , which with γ 1 = 0 yields γ r = r − 1.If we consider equation (33) then the summation is dominated by the first term and so p 1 = p 0 − az 2  1 .Finally, the leading order terms of equation (31) yield which, with the leading order terms from equation (42), implies These results have been confirmed by numerics [24], measurements of the gradient of a log-log plot of x 3 and x 4 against time yield 1.99 and 2.97 respectively.Thus as t −→ ∞ we have x r (t) ∼ (2ax 0 ) r/2 /2(2p 0 − x 0 ) r−1 t r−1 for r > 1. (47) The borderline case, 2p 0 = x 0 , is such that there is just enough inhibitor to poison all the clusters formed in the system, but it is not obvious how fast this will occur.Numerical simulations suggest that the following scalings are appropriate for the problem with ξ = 1/2, and this will be confirmed by the following analysis.We assume no a priori knowledge of χ 1 , χ r , z r , z 1 or ξ, only ξ > 0. The leading order terms in (31) give x 0 − 2ax 2 1 = 0 hence z 1 = x 0 /2a.Substituting equations ( 48)-( 50) into (33) yields which to leading order gives p 0 = kp 1 z 2 .Balancing the next order terms gives ξ = 1/2 and so 1 2 With the scaling parameter we now consider (32) The leading order terms occur at O t −(r−2)/2 yielding a recurrence relation with solution Keeping terms to O(t −(r−1)/2 ) gives a recurrence relation for χ r , that is where, for ease of analysis, we define the constant Solving the recurrence relation yields Finally, we consider (31), wherein the O(1/ √ t) terms yield We now wish to solve equations ( 52), ( 56) and (58) to find p 1 and χ 1 .Equation ( 52) gives which, when combined with equation ( 56) and (58), yields where the Υ term arising from p 2 on the LHS of equation ( 59) cancels exactly with that arising from the p 1 χ 2 /z 2 term on the RHS; thus there is no need to calculate the value of Υ explicitly.Additionally we have and so the leading order terms have been resolved, however, p 2 remains undetermined since Υ has not been evaluated (that calculation requires higher order equations and is omitted here).In Figure 2A we plot the numerically-determined poison concentration over time; it shows an indefinite increase and using a log-log scale we have confirmed that the timedependence scales with t 1/2 .If we assume that p = p rmnum 1 t 1/2 then by considering the last data point (t = 10 9 ) we obtain p num 1 = 1.9167 which is in good agreement with the analytical result of p 1 = 1.9168 from (60).Additionally we have numerically confirmed that the second-order term decays according to t −1/2 as predicted by the analysis.The monomer concentration is predicted to tend to a finite limit, and this is shown in Figure 2B.The analysis has predicted this limit to be x eq 1 = 0.6124 and this agrees with the numerical simulation.Finally we compared the numerical value of χ num 1 = 0.2477 with the analytical prediction from (61), χ 1 = 0.2477, and they agree.In this case monomers enter twice as fast as poison and while the monomer concentration saturates to x 1 −→ ax 0 /2, the poison accumulates according to p ∼ (ax 0 /2) 3/4 t 1/2 / √ ak, thus the clusters x r for r > 1 are scarce, and these concentrations decay algebraically, the most numerous being the dimers which decay according to x 2 ∼ (x 0 /2a) 1/4 / √ kt.
3.4.General a r , b r , k r , with x 0 < 2p 0 Previous analysis has focused on the special case of size-independent rates (a r = a, b r = b and k r = k), which has been solved for all parameter regimes.In this Section we solve the case of x 0 < 2p 0 but for general rates, that is, size-dependent non-zero a r , b r and k r .The following analysis is similar to that of Section 3.2.One steady-state is x r = 0 (r > 1) with x 1 tending to a finite value, while the poison concentration grows linearly.We thus assume x r ≈ z r /t γr and p = p 1 t; further assuming that γ 1 = 0 in the final timescale.The ẋr equation yields and thus we deduce the same balance as before so that γ r = r − 1.The leading order balance for the above equation will be the first and last term on the RHS and so To calculate p 1 we consider the leading order terms in the equation for ṗ (33) and obtain, as before, p 1 = p 0 − a 1 z 2 1 .We calculate z 1 by balancing the leading order terms in equation ( 8), thus z 1 = x 0 /2a 1 , and so p 1 = p 0 − x 0 /2, and so we can expect this solution to fail if 2p 0 ≤ x 0 .Thus as t −→ ∞ we have the solution p ∼ (p 0 − x 0 /2)t and . This result is independent of the choice of rate coefficients and poisoning susceptibility, assuming that these rates are non-zero.

Summary
In this Section we have considered a system with poisoning and monomer influx, but as yet without competition for monomers, that is, we have allowed only one morphology of cluster to form.We find three different cases, which depend on the relative levels of monomer to poison influx.Firstly, in Case A the poison influx is relatively small; in this case the concentration profile tends to a steady-state solution and the poison level to a finite value and these have been calculated explicitly.Numerical simulations reveal that the concentrations of clusters and of poison decay exponentially to their final values.We additionally note that the average cluster size is unaffected by the susceptibility of poisoning (k), but the final level of poison does depend on k.Secondly, in Case B the poison influx is relatively high compared with the monomer influx.This results in linear, unlimited, growth of the poison concentration (p(t) ∼ (p 0 − x 0 /2)t) while all the cluster concentrations tend to zero in the large time limit, with the exception of the monomer concentration which tends to x 1 = x 0 /2a.The timescale over which the concentrations tend to zero is found to vary according to x r (t) ∼ O(t −(r−1) ) for r > 1.This can be generalised to size-dependent aggregation, fragmentation and poison susceptibility rates where, in general, x 1 −→ x 0 /2a 1 , x r −→ 0 and p ∼ (p 0 − x 0 /2)t as t −→ ∞.In Case C, we consider the case where the monomer and poison influxes are balanced, the poison grows without limit, according to p ∼ (ax 0 /2) 3/4 t 1/2 / √ ak, the monomer concentration again tends to x 1 = x 0 /2a and the larger cluster concentrations tend to zero in the large time limit.However, the larger cluster sizes tend to zero at a slower rate than in Case B, with x r ∼ t (r−1)/2 .
The intuitive reason for a critical point existing at 2p 0 = x 0 is that since we do not allow monomers to be poisoned, dimers must form before the cluster can be poisoned, so with 2p 0 < x 0 the influx of poison is insufficient to poison all the monomers flowing into the system.However, with 2p 0 > x 0 then there is sufficient poison influx to poison all the mass as dimers and maintain a build up of poison.When 2p 0 = x 0 , most of the monomers added to the system forms dimers and is poisoned by the inhibitor, however a few trimers are formed (these have a concentration which is O(1/ √ t) smaller than the dimer concentration).The poisoning of a trimer leaves a slight excess of inhibitor, thus the concentration of inhibitor rises at a rate proportional to O(1/ √ t), and so p ∼ O( √ t).

Competition
We now allow competition between two morphologies of cluster growing from monomers, while the inhibitor only affects one morphology, say the x r clusters, as given by ( 9)- (13).There are four types of behaviour observed depending on the aggregation and fragmentation coefficients and the relative size of 2x 0 to p 0 ; these are discussed after a brief description of the case p 0 = 0, that is the case of competition without poisoning.

4.1.
No poisoning, p 0 = 0 In the absence of poison the input of monomers must be balanced by either the x r or y r clusters tending to the borderline equilibrium solution, x r = x 1 or y r = x 1 ; unless a/b = α/β in which case they will both tend to this state.If a/b > α/β then that is, the x r -clusters are the more thermodynamically stable and will tend to the borderline equilibrium solution, and the y r -clusters will tend to an equilibrium solution which decays as r → ∞.If a/b < α/β then this situation is reversed so that Here x r -clusters are more thermodynamically stable than the y r -clusters and the input of monomers dominates that of poison.In the parameter regime a/b > α/β and without poison we expect the x r -clusters to approach the borderline equilibrium solution x eq r = x 1 and the y r -clusters to form an equilibrium solution y r = θr−1 y x 1 with θy = αx 1 /β < 1. Figure 3 shows the large-time solution of the system both with and without poison.The addition of poison alters the structure of the x r concentrations and the y r concentrations remain in equilibrium, though with x 1 and θy reduced; from x r = x 1 to x r = θr−1 x x 1 with θx < 1 as in Case A (see Section 3.1, equations (39), ( 40) and ( 41)).In Figure 3 we also plot the large-time solution for the system with an influx of poison.By poisoning the x r clusters sufficiently, the y r become the more stable species as expected, although we note that due to the poisoning of the x r clusters the concentrations of the y r clusters has also fallen; this is due to the reduced monomer concentration, x 1 .In general, there will be a critical value of inhibitor influx, p crit 0 , such that θx = θy so that x r = y r ∀r; by combining θy = αx 1 /β and equations ( 39) and (40) we deduce that p crit 0 satisfies βb α which clearly has a unique root for p crit 0 between 0 and x 0 /2.Here the x r -clusters remain more stable than the y r -clusters, but there is now an abundance of inhibitor.Case II differs from Case I only by the increase of inhibitor in the system, thus without poison the concentration profiles will be identical to those shown in Figure 3.With an increased amount of poison, p 0 > x 0 /2, the y r -clusters remain in equilibrium (y r = x 1 (αx 1 /β) r−1 ) while the x r -clusters behave as described in Case B (Section 3.2), thus x r −→ 0 for r > 1 and x 1 −→ x 0 /2a as t −→ ∞.In this case the poison concentration rises linearly without limit, p ∼ p 0 t, and the x r are completely inhibited from growing, and so the y r clusters remain as the only stable form of cluster in the system, as shown in Figure 4.In summary, as t −→ ∞, we have   In this case the y r -clusters are more stable than the x r -clusters and there is an abundance of monomer over inhibitor.In Figure 5 we plot numerical results of the case with no input of inhibitor, showing that the y r -concentrations approach the borderline equilibrium solution y r = x 1 = β/α.From Section 2 we recall that the monomer concentration adjusts such that αx 1 /β = 1 in the large time limit.The x r concentrations form the equilibrium solution x r = (ax 1 /b) r−1 x 1 = (aβ/bα) r−1 β/α.The addition of any amount of poison causes the monomer concentration to fall and thus the y r -clusters to revert to an equilibrium solution which decays with increasing r, given by y r = θr−1 y x 1 with θy < 1.The condition 2p 0 < x 0 implies that the x r -clusters adjust to a steady-state solution given by x r = θr−1 x x 1 .Since the y r -clusters are in equilibrium, the analysis of Case A (Section 3.1) holds exactly and θx , θ and p final are given by equations (39), ( 40) and (41) respectively, as shown in Figure 5.This is a surprising result.Without poison the y r clusters form increasingly large clusters and are by far the more stable configuration.While it is true that with poison the y r clusters are still the more stable morphology, the form of the solution for the y r clusters is altered considerably, from the borderline equilibrium in which the mass of y r -clusters grows without bound ( r ry r = O(t)) to an equilibrium solution with a finite mass; a change induced by poisoning the less stable cluster type (x r ).4.5.Case IV: a/b < α/β, x 0 < 2p 0 In the final case, y r -clusters are more stable than x r -clusters and the influx of poison dominates that of monomers.If no poison is present then we still expect the long time results shown in Figure 5, however, we now investigate the effect of a large amount of poison on the system, that is 2p 0 > x 0 .The y r -clusters form an equilibrium solution y r = θr−1 y x 1 and so the analysis of the x r clusters follows exactly that of Case B (Section 3.2); thus in the large time limit we expect x r −→ 0 as t −→ ∞ for all r > 1 and the monomer concentration to be given by equation (43).In Figure 6 we plot the numerical x r and y r concentrations in the large time limit and note that equation (43) predicts a monomer concentration of x 1 = 0.6455 which agrees well with the numerical simulation.Again it is interesting that poisoning the less stable clusters causes a dramatic change in the form of solution of the more stable cluster type.

Conclusion
In this paper the Becker-Döring model has been generalised several times to build up a complicated model of competition for an influx of monomers, between two types of morphologies, with one morphology being susceptible to an inhibitor which is also being introduced at a constant rate.Numerical simulations have been performed to confirm the analytical results obtained, after truncating the model at a large finite size r = N , typically 500.Initially we considered the Becker-Döring model with a constant influx of monomers.At small times the clusters form an equilibrium-like structure with θ = ax 1 /b < 1.We found that as t −→ ∞, x 1 −→ b/a so that the borderline equilibrium solution x r = x 1 = b/a is approached.We calculated the similarity solution which is valid in the large-time limit (30).
The next generalisation added to the model is allowing a constant influx of inhibitor.When a cluster is poisoned it is effectively removed from the system, since it plays no further role in the kinetics of cluster growth.The solution is found to depend on the relative influx of monomer (x 0 ) to inhibitor (p 0 ), with a special case when 2p 0 = x 0 .In the case 2p 0 < x 0 we found that the poison tends to a finite concentration while the x rclusters tend to a steady-state solution.In contrast, if x 0 < 2p 0 then the concentration of inhibitor increases without limit, linearly in time, and x r −→ 0 for r > 1, with x r = O(t −(r−1) ) as t −→ ∞.In the borderline case 2p 0 = x 0 , the system again tends to the solution x r −→ 0 for r > 1 as t −→ ∞ and the inhibitor concentration tends to infinity as t −→ ∞.In this case the appropriate scalings are p ∼ t 1/2 , x 1 ∼ x 1 and x r ∼ 1/t (r−1)/2 ; and we calculated the leading order terms in all cases.In the case x 0 < 2p 0 the scalings hold for general a r , b r and k r , calculating the final solution and the timescale over which this is approached.If clusters of size r < m were immune from poisoning, then we expect a critical input ratio at x 0 = mp 0 , with all cluster concentrations decaying to zero if x 0 < mp 0 and a non-zero steady-state if x 0 > mp 0 .
The final generalisation to be added to the model is to allow competition for monomers, that is, where monomers aggregate into one of two types of cluster, x r or y r , but we allow only the x r -clusters to be poisoned.In the absence of a poison influx, we find that the less thermodynamically cluster forms an equilibrium solution which decays with increasing r, while the more stable morphology approaches the degenerate equilibrium solution in which the concentrations of all sizes are equal x 1 .In the presence of inhibitor, we first assume that the y r -clusters are the less stable, then poisoning the x r -clusters results in a steady-state solution for the x r -clusters if 2p 0 < x 0 , and if x 0 < 2p 0 then x r −→ 0 for all r > 1 (identical results as previously obtained without competition).However, if, without poison, the y r -clusters are more stable and form the degenerate equilibrium solution y r = x 1 , then poisoning the less stable cluster, (x r ), results in the y r -clusters forming an equilibrium solution which decreases at large cluster sizes.The size distribution of the more stable morphology is significantly altered by poisoning the less stable cluster; this is not an intuitive result.A naive interpretation of such an observation would lead one to assume that the inhibitor was acting primarily on y r -clusters, or possibly on both the x r -and the y r -clusters.In the unpoisoned case the less stable morphology (x r ) reaches its equilibrium and then all input monomers form the y-morphology, whereas with an inhibitor a steady-state is reached in which all input monomers form the x-morphology and are then poisoned.
These results have important applications in polymorph prediction, and in explaining why the predicted most stable morphology of a crystal is not always the one observed in nature.In cases where two polymorphs have a similar thermodynamic stability, the presence of other chemical species may have a stronger inhibiting effect on the growth of the more stable morphology, thus allowing the less stable morphology to dominate.We hope to generalise this work to the more complex cases where the aggregation and fragmentation rates are size-dependent.This occurs, for example, in classical nucleation theory (see Lewis [25] for details), where statistical mechanical models of the growing crystal nucleus are formulated in terms of a surface energy and a bulk energy.As well as systems which can form multiple morphologies of crystal, models such as these could also be applied to the problems of protein crystallisation where amorphous solids as well as crystals are often produced, see Kam [26] and Weber [27] for further details.

Figure 1 .
Figure 1.A) Difference between the numerical value of p(t) and the analyticallydetermined p final for times 0 < t < 800; illustrating exponential convergence.B) Comparison of the numerical solution x r (2000) and the equilibrium x r = θr−1 x bθ/a.Parameters used were x 0 = 1, p 0 = 0.4, a = 0.5, b = 1 and k = 0.25.

Figure 2 .
Figure 2. The special case x 0 = 2p 0 : (A) the poison concentration against time; (B) difference between the monomer concentration and its predicted limit, against time.The parameters used were a = 2, b = 1, x 0 = 1.5, p = 0.75 and k = 0.25.