A Thermodynamic Model of Microtubule Assembly and Disassembly

Microtubules are self-assembling polymers whose dynamics are essential for the normal function of cellular processes including chromosome separation and cytokinesis. Therefore understanding what factors effect microtubule growth is fundamental to our understanding of the control of microtubule based processes. An important factor that determines the status of a microtubule, whether it is growing or shrinking, is the length of the GTP tubulin microtubule cap. Here, we derive a Monte Carlo model of the assembly and disassembly of microtubules. We use thermodynamic laws to reduce the number of parameters of our model and, in particular, we take into account the contribution of water to the entropy of the system. We fit all parameters of the model from published experimental data using the GTP tubulin dimer attachment rate and the lateral and longitudinal binding energies of GTP and GDP tubulin dimers at both ends. Also we calculate and incorporate the GTP hydrolysis rate. We have applied our model and can mimic published experimental data, which formerly suggested a single layer GTP tubulin dimer microtubule cap, to show that these data demonstrate that the GTP cap can fluctuate and can be several microns long.


Introduction
Microtubules are dynamic filaments that perform essential functions in eukaryotic cells including nuclear and cell division and intracellular transport. A microtubule is a cylindrical assembly of tubulin dimers which are composed of a and b-tubulin subunits. These dimers associate head to tail to form protofilaments and usually 13 protofilaments associate laterally to form the wall of the microtubule. The protofilaments are slightly shifted with respect to each other and after one full turn of the microtubule there is a total shift of 1.5 dimers at the seam, the join of the first and thirteenth protofilament in the microtubule [1].
The tubulin dimer has two GTP binding sites a nonexchangeable GTP site at the interface of the a and b subunits and an exchangeable site between the dimers in a protofilament. At the exchangeable site GTP hydrolysis occurs 250 times faster when the GTP tubulin dimer is bound within the microtubule compared to when it is in free solution [2]. GTP-tubulin dimers have a slightly altered conformation compared to GDP-tubulin dimers and as a result only the former can assemble into microtubules. When hydrolysis of GTP takes place within the microtubule it is the neighbouring dimer interactions that hold the GDP dimers in place [3].
The asymmetry of the tubulin dimer is translated into the microtubule as it assembles and the exposed b tubulin end is called the plus end and the other end the minus end. Each of the plus and minus ends have different properties in respect of structure and growth [4]. When sufficient free GTP tubulin dimer is present microtubules will grow and the result is a cap of GTP tubulin dimer on the microtubule. As the growth rate slows the cap is lost as the GTP is hydrolysed to GDP in the tubulin dimer. As the binding energy of GDP-tubulin dimer is lower than that of GTPtubulin the microtubule undergoes rapid shortening. This selfassembly and disassembly of microtubules is known as dynamic instability [5].
Modelling microtubule dynamics is giving further insight into the manner of microtubule assembly and disassembly. So far the best approach to model microtubule dynamics has been to use Monte Carlo simulations which were first performed by Chen and Hill [6,7] initially for a single protofilament then progressing to a 13 protofilament microtubule. This model had 17 parameters some of which were chosen so as to reproduce the experimental values of Mitchison and Kirschner [5]. A few years later, Bayley and colleagues [8] proposed a similar model, except that GTP tubulin dimers were assumed to hydrolyse spontaneously once embedded inside the microtubule. Such models lead to a GTP tubulin cap that is only one heterodimer long. However, Van Buren et al [9] have introduced a model where the number of parameters was reduced to 4 parameters per microtubule end by relating the tubulin attachment and detachment rates to the binding energy of tubulin heterodimers.
In this paper we derive a thermodynamic model for microtubule dynamics and use this model to perform Monte Carlo simulation where we include the contribution of water to the entropy of the system and we fit all the parameters of our model to published experimental data including the hydrolysis rate of GTP. We consider both the + and the 2 ends. Moreover, using our model we have reinterpreted existing experimental data [10] that were used to predict a short GTP cap to show that the cap can be several microns long and dependent on the concentration of free GTP-tubulin dimer, which is now consistent with recent observations [11].

Theoretical model: thermodynamics
In this paper, we consider the microtubule as a lattice with a 1.5 dimer shift at the seam (Fig. 1). We also view the cap as having two components: a crown consisting of incomplete protofilaments and the core that forms the body of the complete microtubule and includes GTP tubulin dimers (Fig. 2).
The polymerisation of microtubules involves 3 types of reactions. First of all, GTP-tubulin dimers, which we denote by T GTP , can attach or detach from the tip of the microtubule. GDPtubulin dimers, T GDP , on the other hand, can only detach from the microtubule tip. Then, T GTP inside the microtubule can hydrolyse into T GDP .
To describe the polymerisation of microtubules one must consider all the possible configurations that a microtubule can take. Each microtubule is characterised by the length of each protofilament as well as the type of every dimer in each of these protofilaments. The number of possible microtubule configurations, which we call Q, is thus extremely large and instead of using a sophisticated labeling system to describe each of these configurations, we have decided to label them formally using a single index J that takes Q different values. As we will not need to consider the details of each configuration, this formal parametrisation has the advantage of being both simple and sufficient for what we want to do.
We then denote by M J a microtubule in the configuration J. Each time a tubulin dimer detaches from or attaches itself to a microtubule, or each time a GTP-tubulin dimer hydrolyses to a GDP-tubulin, the configuration of the microtubule changes.
The polymerisation of microtubules involves 3 types of reactions. First of all, GTP-tubulin dimers, which we denote by T GTP , can attach or detach from the tip of the microtubule. GDPtubulin dimers, T GDP , on the other hand, can only detach from the microtubule tip. Then, T GTP inside the microtubule can hydrolyse into T GDP . These three reactions, which all have their own reaction rates described later, can be summarised as follows: where Pi is an inorganic phosphate. We have added here, for completeness, the self regulating process transforming GDPtubulin into GTP-tubulin. Note also that we have enclosed within parentheses the back reaction arrows for the last 3 equations in (1) because, while their reaction rates are very small, they must still be considered for a full thermodynamic description of the system. In each of these equations, the indices J, J', J'' and J''' are related to each other but we still have a very large number of simultaneous reactions involving different concentrations of various microtubule configurations. Equation (1) is a formal system of equation corresponding to a very large number of chemical equations, one for every combination of related indices. The last equation in (1) symbolises the self regulatory process that controls the GTP-tubulin concentration in the cell. For ie in vitro experiments, the number of microtubules is usually small enough that the GTP-tubulin concentration remains constant. For both types of experiments we will thus assume, in what follows, that the GTP-tubulin concentration is constant.
It is important to realise that each microtubule can transform into many other microtubule configurations and that each microtubule can also be obtained from many other microtubule configurations. A specific microtubule configuration M J can occur on both sides of the first 3 reactions in (1) for a small set of explicit values of J, J', J'' and J'''. Equilibrium for (1) can then be achieved because the microtubules decaying, say by hydrolysis, are recreated through a combinations of all 4 equations.
We now consider a volume V containing N J moles of microtubules in configuration J, N t moles of free GTP-tubulin dimers, N w moles of water, N d moles of free GDP-tubulin dimers and N p moles of Pi. A priory, we thus have Q different types of microtubules, each with their own concentration N J . Moreover, we know that in physiological conditions, N J =N w , N t =N w , N d =N w and N p N w are always very small.  The Gibbs' energy of the microtubule solutions is given by [12] G(N J1 ,:: where u T ð Þ, s T ð Þ and v T ð Þ are, respectively, the energy, the entropy and the volume of one mole of the different constituents and each of these quantities depends on the temperature T. As a matter of fact in (2), for each solute, i.e. not for water, we have s X T ð Þ~s s X T ð ÞzR, wheres s X T ð Þ is the entropy of the corresponding solute while R is the entropy gained by the solvent, water, as the solute is added to the solution. s X T ð Þ ands s X T ð Þ are thus two valid expressions of entropy but which have different physical origins. The use of one rather than the other will not affect our results as such but it will modify the interpretation of some of the quantities discussed later.
For a perfect solution, the energies u i and the entropiess s i are, respectively, the internal energies and internal entropies of each substance.
In (2), R~N a k where N a is Avogadro's number and k is the Boltzmann constant.
During the association of one tubulin dimer to a microtubule M J to produce a microtubule M J' , as in the first reaction in (1), the Gibbs' energy changes and this variation can be easily computed as Note that the variations of Gibbs' energy during the dissociation, as in the first two reactions in (1), are also given by (3), after changing the signs on the right hand side and, for the detachment of GDP-tubulin, exchanging N t and N d .
We have decided to compute all the quantities normalised to single molecules rather than 1 mole because we want to derive a model based on a Monte Carlo simulation of the chemical reaction which deals with only one molecule at a time.
Using the fact that log N J {1=N a ð Þ &log N J we have Note that an equation like (4) can also be used, up to a sign, for the detachment of any type of tubulin dimer. A similar expression can be used to compute the variation of the Gibbs' energy during the hydrolysis, but as we will not need it we do not present it here.
Denoting, respectively, the energy, entropy and volume per molecule as U~u=N, S~s=N and V~v=N, we define the variation of energy, entropy and volume during the reaction as: Notice that the variation of energy DU J,J' is the energy necessary to detach a single tubulin dimer from the microtubule in configuration J to obtain one in configuration J' in the considered solution. In an ideal solution, this is just the binding energy of the molecules while in a non-ideal case, DU J,J' will also include the variation of all electrostatic/chemical interactions and of the average configuration changes of all the tubulin dimers, bound or not, and of the solution itself. This could, in theory, be computed from the full Hamiltonian of the system which includes all the interactions ( ie chemical, or electrostatic) between the different substances present in the solution. A given u i in (2) will thus be the sum of the internal energy and of half the interaction energy with the other constituents. By symmetry, each interaction term will occur twice, hence the need to half them. One should also take into account the presence of other solutes, such as ions or proteins , always present in experimental assays or in in vivo solutions. Such solutes do not play a direct role in the dynamics of the microtubule, but they can affect indirectly the binding energies of the tubulin dimers and hence, as we will see, have an indirect impact on the microtubule dynamics.
One could, ideally, evaluate these variations of energies by solving the Schrödinger equation describing the system, but in practice this is impossible given the number of electrons involved. An alternative would be to use semi-classical methods to approximate the solutions of the Schrödinger equation. Then, the interaction between molecules will be described by a collection of tailored chemical and electrostatic potentials describing the interactions between the tubulin dimers, water molecules and the other substances present in the solution [13].
In practice, the cytoplasmic solution of microtubule and tubulin is not an ideal solution in the thermodynamic sense as there are non zero forces between all the molecules involved. DU J,J' thus contain contributions which, in principle, could be computed ab initio, but which are usually described by the activity parameters for the reaction. In practice we will fit the values of DU J,J' to experimental data without attempting to derive them.
Note also that

Microtubule at Equilibrium
When microtubules are studied experimentally either in ie in vitro or in vivo experiments, their dynamics is usually not at thermal equilibrium. Thermal equilibrium is reached when a solution has a very large number of microtubules polymerising and depolymerising according to (1) and where the concentration of all the solutes remain, on average, constant. As we will see below, this critical concentration is very special and it will allow us to link together the Gibbs' energy, computed above, to the attachment and detachment rates of tubulin dimers and to the value of the critical concentration.
The thermodynamic equilibrium is characterised by the condition DG J,J'~0 and the chemical potential Dm J,J' is defined as To evaluate this variation of the chemical potential, we assume that the entropy of a microtubule composed of n J tubulin dimer is n J times the entropy of a single tubulin dimer:s s J~nJs s t . This assumes that the tubulin dimers have some vibrational degrees of freedom when they are attached to a microtubule. As a direct consequence, the entropy does not change during the binding process and DS S J,J' &0. This in turn implies that, for a perfect solution DS J,J' &{k. If our approximation were incorrect, then DS S J,J' v0. DS S J,J' would also be non-zero for a non perfect solution, but for the sake of simplicity, we will take DS J,J'~{ k at this stage and we will show later that, if DS J,J' were different, the model, as it stands, would not be affected. Only the derived value of the longitudinal binding energy of a tubulin dimer would have to be changed.
In a liquid solution we have v J' &v J zv t which implies that the volume of the system does not change significantly during the binding process either and DV J,J' &0. In our model, so far we have not taken into account the fact that the entropy of the system might change as a result of electrostatic interactions or because of some geometrical effects. For example, when a tubulin dimer detaches itself from the microtubule it exposes a larger area on which water molecules can cluster themselves on the protein. We would expect this to decrease the entropy.
We will assume that such entropy effects are small, but if they were not, they could be added to DS J,J' which would then be non zero. Similarly, because of the electrostatic interactions between the dimers and the solution, DV J,J' might not be small enough to be negligible. In either case, at constant temperature and pressure, non zero factors for TDS J,J' and PDV J,J' would just be constant contributions to m J,J' which would not affect the modeling results but, as we shall see later, would alter the interpretation of the binding energy.
Before we proceed further, we rewrite (7) so that it expresses the relative number of moles of the different microtubule types as well as GTP-tubulin at equilibrium as a function of the tubulin dimers binding energies: We should stress that this relation is only valid at the critical concentration, which corresponds to the thermodynamical equilibrium of the system, and where, moreover, the hydrolysis rate is small compared to the polymerisation and depolymerisation rates. Note also that (8) is dimensionally balanced.

Microtubule (De)-Polymerisation Rates
Considering the same volume V as above, we must now analyse the rate at which the microtubule (de-)polymerises. First of all, we define k T {;J',J and k D {;J',J as the dissociation rates per microtubule of, respectively, a GTP-tubulin or GDP-tubulin dimer for a microtubule going from state J' to J. k z is then defined as the first order rate at which a microtubule polymerises to another state at a given concentration of free GTP-tubulin. GDP-tubulin does not polymerise, so k z for GDP-tubulin is zero. As the mean free path of a tubulin dimer in a water solution at body temperature is very short when compared to its size, we can assume that k z for GTPtubulin does not depend on the the microtubule state. In other words, we assume that there are no geometrical factors such as there would be if we considered a gas. We also define r h as the hydrolysis rate of GTP-tubulin into GDP-tubulin inside the microtubule and n pf as the number of proto-filaments.
The rate of change of N J is given as the sum of all the transitions creating a microtubule in state J less the sum of all the rates of transitions from which state Jcan decay. For the sake of convenience, we define the following sets of microtubule states: N I T z : the set of n pf states which can be obtained from state J by adding a GTP-tubulin dimer. N I D z : the set of n pf states which can be obtained from state J by adding a GDP-tubulin dimer. N I T { : the set of states which can be transformed into the state J by removing one GTP-tubulin dimer. N I D { : the set of states which can be transformed into the state J by removing one GDP-tubulin dimer. N I hz : the set of states from which the state J can be obtained through the hydrolysis of one GTP-tubulin dimer. N I h{ : the set of states into which the state J can transform through the hydrolysis of one GTP-tubulin dimer.
Note that together, the two sets I T { and I D { have a total of n pf different states.
Using these definitions, the variation of N J , the number of microtubules in configuration J, is given by depolarisation or T GTP hydrolysis from microtubules in other configurations, reduced by the number of microtubules M J that are converted to other microtubule configurations also by polymerisation, depolarisation or hydrolysis. At thermodynamic equilibrium, dNJ dt~0 and, following [9], we can use that special configuration to determine a relation between k z , k T {;Jz,J and k D {;Jz,J . In the case of a pure GTP-tubulin microtubule when r h~0 , it is straightforward to show that is a solution of (9), where we have introduced the parameters k z;eq , N J;eq and k T {;Jz,J;eq defined as the values of k z , N J and k T {;Jz,J , respectively, at thermodynamic equilibrium. Eq. (10) states that if no hydrolysis takes place then microtubules are made entirely out of GTP-tubulin and so the polymerisation rates for J?J' and depolymerisation rates for J'?J are the same.
When the hydrolysis rate, r h , has a value comparable to the other dynamical rates, then finding a solution to (9) is very hard, and one cannot derive any expression like (10). This is because the ratios of k z and k { for different J and J' are not directly related to their relative concentrations anymore, but the rates and concentrations for different configurations are all interdependent.
If the hydrolysis rate is non-zero but small, (9) has several residual terms, but we can assume that they are all very small. Indeed, if r h is small compared to k z and k { , the GTP cap will be large, and there will be many states from which and into which the state J can hydrolyse. Moreover, on average, the number of these states will be very similar and they will have similar probabilities to be generated. We thus see that the 2 terms in the last sum in (9), proportional to r h , will mostly cancel each other out.
If the hydrolysis rate is small compared to k z , there will be very few GDP-tubulin dimers at the tip of the microtubule. This implies that there will be few microtubules in the configurations I D z and I D z and that the sums over these configurations in (9) will be small. This will indeed be confirmed later by our simulations. In this case (9) is nearly correct and we can use it as a good approximation to the relation between k z and k T {;Jz,J for the purpose of our study.
Note that in Eq. 9 k z;eq , k T {;eq , k D {;eq and r h are all rates expressed in units of ''per second per dimer'' and they can thus all be compared with each other. To be able to use Eq. 10 we must have r h vk z;eq , k T {;eq and k D {;eq . (k z;eq is also expressed in units of per second per dimer because it is defined as the attachment rate at the critical concentration of free GTP-tubulin).
To estimate k T { as a function of k z , we ignore the hydrolysis rate and consider it as a small correction in the Monte Carlo simulation. This approximation will be justified if we show that, indeed, r h is small.
Within that approximation, we can use (8) and (10) which establishes a relation between the dissociation and association rates of microtubules and free tubulin at the thermodynamical equilibrium. As shown in (11) the ratio is related to the difference in energy of the 2 microtubule configurations and to their relative concentration at the critical concentration.
The dissociation rate of any dimer close to the tip of the crown can be computed using the expression k z;eq k T {;J',J;eq~r where r w~5 5:56 molar is the molar concentration of water and r t;eq is the molar critical concentration of GTP-tubulin. We can thus use (12) to determine k T {;J',J;eq as a function of k z;eq and of the binding energy difference which depends on J and J'.
It is also important to stress that the derivation of (12) is only valid at the equilibrium, i.e. only at the critical concentration r t;eq . Moreover (12) assumes that r h is small compared to k z;eq and k T {;J',J;eq . When the concentration r t is different from r t;eq , we do not have statistical equilibrium, but the only parameter that changes is k z and we can write where d is the relative tubulin concentration defined as the tubulin concentration normalised to the critical concentration: To perform a Monte Carlo simulation outside the equilibrium conditions, we use the fact r h does not depend on r t and that the dissociation parameters k T {;J',J~k T {;J',J;eq are the same as at the equilibrium and can thus be computed using (12). This is equivalent to the assumption that the detachment rate does not depend on the free tubulin concentration.
We must stress here once more that eq. (12), which relates the depolymerisation rates of the microtubules to their variation of internal energy, is only valid if the hydrolysis rate r h vanishes, or, as an approximation, if it is small compared to the first reaction at equilibrium. If r h is too large, then there is no simple expression like (12). Instead, the depolymerisation rates and binding energies of all the possible states would depend on each other. Fortunately this turns out not to be the case.

Monte Carlo Simulation
Our model involves a Monte Carlo approach to simulate the attachment and detachment of tubulin dimers to and from a microtubule. At any given time, several events can take place with different probabilities: N A GTP-tubulin dimer can attach itself at the tip of the microtubule at the rate k z given by Eq. 13 . N A GTP-tubulin dimer can hydrolyse into GDP-tubulin at the rate r h . For the + end, as the exchangeable GTP of the GTP-tubulin dimers at the tip of any protofilament are not embedded inside the microtubule they hydrolyse very slowly and so we can put r h~0 for these. All the other dimers hydrolyse at the same rate. N Any dimer, GTP-tubulin or GDP-tubulin, can detach itself from the microtubule. The detachment rates k T { and k D { are obtained from Eq. 11 , k {~kz;eq r w r t;eq e Dm J,J' kT , where Dm J,J' is the variation of chemical potential during the detachment, which depends on the microtubule configuration and on the position and the type of dimer. N Strips of dimers are also allowed to detach together. In that case the variation of chemical potential Dm J,J' correspond to the total binding energy of the attached dimers. In practice, as jDm J,J' j tends to be large, this occurs relatively rarely, except for dimers with no neighbours on one side.
As the detachment rate k T {;J',J is independent of d it can be evaluated using Eq. 11 . To do this, we must determine Dm J',j . Note that to use Eq. 11 we must also know r t;eq which has been determined experimentally by Walker et al. [14].
To evaluate Dm J,J' we use Eq. 7 and we observe that the binding energy of tubulin dimers can be split into the lateral and longitudinal binding energies which we denote, respectively, e l and e L in units of kT. To follow the conventions of Van Buren et al [9], we define G L~eL {1, because DS J,J'~{ k, and G l~2 e l (see Fig. 1). Then we can write where the superscript T and D refer, respectively, to GTP and GTP-tubulin and where Dn J,J' is the number of lateral monomer bonds that must be broken to go from configuration J' to configuration J. Note that by taking into account the entropy of water in Eq. 8 , the binding energy e L is 4kT larger than the binding energy derived in [9] for an effective gas. Moreover, if we decide to modify the assumptions that we have made regarding the entropyS S and assume DS S J,J'~k s where s is a quantity that one would have to determine, then we would have G L~eL {1zs. As G L is the parameter that is fitted from experimental data, we see that our model is not affected by our assumptions, but that the value of the longitudinal binding energy e L , as a side product, is. When several dimers detach themselves together from the tip of a protofilament, Dm D J,J' and Dm D J,J' are still given by Eq. 15 and Dn J,J' can be larger than 4. In practice, however, the detachment of several dimers when Dn J,J' is large is a very rare event as such events are exponentially suppressed.
To evaluate k D {;J',J we also use Eq. 11, taking the same value for G L and k z;eq and for G l we take the value that fits the GDPtubulin depolymerisation rate which we call G D l . In doing so we exploit the fact that GDP-tubulin and GTP-tubulin are very similar dimers and we assume that the difference of binding energy comes from the lateral bounds; as GDP-tubulin is curved [15], it is less bound than GTP-tubulin.
To evaluate k T {;J',J and k D {;J',J we must thus determine the three binding parameters G L , G T l and G D l . Ideally, one would like to derive them by performing an ab initio computation of the binding energies between tubulin dimers but this would be extremely difficult and, instead, we have determined their values using the microtubule dynamics experimental data [14].
As stated above, the rates k T {;J',J , k D {;J',J , k z and r h for the four types of events that takes place in the model are all expressed in the same units and can thus be used in a first reaction Monte Carlo simulations [16,17] as follows. For any event with rate k, we computed the time t~{log p ð Þ=k where p is a random number in the range 0,1½ and we do this for the following events: N For every GTP-dimer in the microtubule, we consider its hydrolysis at rate r h and its detachment together will all dimers attached longitudinally above it with rate k T {~k z;eq r w r t;eq e Dm T J,J' kT . N For every GDP-dimer in the microtubule, we consider its detachment and its hydrolysis together will all dimers attached longitudinally above it with rate k T {~k z;eq r w r t;eq e Dm D J,J' kT . N For every tubulin dimer with no other dimer about it, i.e. tubulin dimers at the tip of their protofilament, we also consider the addition of a dimer to the same protofilament with rate k z .
From this large number of random times, we picked the shortest one, t min , and selected the event to which it corresponds. We then implemented that event and increased the simulation time by t min , repeating this procedure for as long as required.
In summary, our model has 9 parameters: the relative free tubulin concentration d, the GTP hydrolysis rate r h and the longitudinal free energy G L , which have the same values for the + and the 2 end as well as the attachment rate k z;eq and the lateral binding energies G T l and G D l which have different values at the + and the 2 end.

Parameter fitting with Experimental Data
To fit the values of the parameters of our model, we have used two sets of experimental data. The first set is made out of the parameters of the microtubule dynamics measured by Walker et al. [14], for both the + and 2 ends, and which are summarised by the expression where R M is the growing rate of the microtubules at relative concentration d, defined by Eq. 14 , while K z;eq and K { are the average microtubule attachment and detachment rates. At the critical concentration, by definition, d~1 and R M~0 and so K z;eq~K{;eq . Walker et al. found that the critical concentration r t;eq~4 :9mM and that, when there were no catastrophes, K z;eq &3:4 dimers/s/protofilament for the + end and r t;eq~5 :3mM and K z;eq &1:8 dimers/s/protofilament for the 2 end. In the catastrophe mode, the depolymerisation rate of GDP-tubulin K D { &56 dimers/s/protofilament for the + end and K D { &70 dimers/s/protofilament for the 2 end. The second set of data comes from the measurement of the time delay before the onset of microtubule depolymerisation performed by Walker et al. [10]. In this experiment, microtubules were polymerised and the growing rates of microtubules measured. Then the solution was washed out by a solution free of tubulin and the time elapsed before the microtubule starts to depolymerise in a catastrophe was measured. The delay observed was roughly equal to 10 seconds, increasing slightly with the initial growing rate.
To ensure that our results are not affected by programming errors, the Monte Carlo simulations for our model were implemented and run totally independently by two of the authors. Their results matched perfectly.
Our first observation is that the growing rate of the microtubule is not linear as a function of the free tubulin concentration but that it is nearly linear in the range covered by the experiments in [14] and [10] (as seen in Fig. 3). This is why we took only 2 values for the concentration to fit the parameters.
To find values of G L and G l to reproduce the experimental data we had to take a k z;eq value larger than K z;eq . To decide which value of k z is best, we had to compare the growing rate of the microtubule as a function of the concentration and compare it to the experimental curve for Eq. 16 in [14].
To determine the best values of the parameters, we have calculated the values of G L , G T l and G D l that reproduce the experimental data of Walker et al. taking r h~0 :1, 0.2 and 0.3 and varying k z;eq between 5 and 9 with an increment of 0.5 for the + end and between 2.5 and 5 with an increment of 0.25 for the 2 end. Note that for a given value of G L , G D l is simply determined by fitting the depolymerisation rate, K D { , of a pure GDP-tubulin microtubule. We have first determined the best values of G L , G D l and G T l for various values of k z;eq and r h by comparing the growing rate and depolymerisation rate values of [14] for d~2 and d~5. In Fig. 4, we present the values of the binding energies G L , G T l and G D l as a function of k z;eq for the + and 2 ends, for r h~0 :3. We have then compared the dynamics and the delay curves computed for these parameters with the experimental data. In particular, the hydrolysis rate r h was determined by matching ( Fig. 5 and 6) the experimental washout delay time of [10]. For this we simulated the dilution process of [10]: a 3 seconds delay followed by 5 seconds during which the free tubulin concentration decreases linearly. We have also chosen G L to be effectively the same for the 2 and the + ends and concluded that the best values    of the required parameters are k z;eq~6 , r h~0 :3, G L~1 2:399, G T l~9 :763, G D l~2 :908 for the + end and k z;eq~3 , r h~0 :3, G L~1 2:566, G T l~1 0:085, G D l~1 :832 for the 2 end. We do not claim to have a good accuracy on these numbers. The error in k z;eq is probably of order 1 and this, by itself, induces errors in G L , G T l , G D l which are relatively large. Unfortunately with the experimental data available, it is not possible to determine G L and G l more accurately (see Fig. 3, 5 and 6).
Note that the individual dimers attachment rates k z;eq is not identical to the average microtubule growing rate K z;eq . The difference comes from the interaction between neighbouring protofilaments which makes the microtubule dynamics non-linear.
Having computed G L and G T l for various values of k z;eq for the + and 2 ends of the microtubule we can read from Fig. 4 which values of k z;eq for the + and 2 end correspond to the same values of G L or G T l . This is shown in Fig. 7 where we see that the two curves are close to each other, indicating that it is G D l that differs between the 2 ends of the microtubule, while G L and G T l are effectively the same for both ends.

Application of the model to microtubule structure
In Fig. 8A and 8B we present two typical snapshots of the plus end of a growing microtubule (from movies S1 and S2 provided in the supplementary information) where the relative density, d, of free tubulin is different but all other parameters i.e. the binding energies and the dimer attachment rates are the same and as derived in the previous section. In Fig. 8A d~2 and in Fig. 8B d~5. Moreover, the microtubules are represented in an unfolded configuration so one can see all 13 protofilaments with the top and the bottom protofilament forming the seam in the folded microtubule. These movies show that the microtubule cap length which is the distance from the tip to the last GTP tubulin dimer is proportional to the concentration of the free tubulin and in these simulations is ca 20nm in movie 1 (d~2) and ca 80 nm in movie 2 (d~5).
In order to analyse the size of the cap we performed a simulation where we sampled the distribution of GTP-tubulin dimer every second. We counted the average number of GTP tubulin dimers from a particular reference point over 100000 samples. This reference point was the base of the microtubule crown i.e. the point at which the 13 protofilaments at the tip are not folded into a complete cylinder (Fig. 2). Moreover, we counted in both directions: into the crown (Fig. 9) and into the core (Fig. 10) of the microtubule cap. Fig. 9 shows the average distribution of GTP-tubulin dimer in the crown for free tubulin concentrations between d~1:25 to d~5 and for one protofilament. As the crown is almost exclusively GTP tubulin dimer, the graphs can be interpreted as representing the time distribution of the protofilament length in the crown and this generally fluctuates from 1 to 7 dimers but can be as high as 10 dimers. For example, we can interpret from the data in Fig. 9 that when d~1:25, the crown is 2 dimers long only 30% of the time whilst for d~5 it is 2 dimers long 50% of the time.
When we measure the GTP-tubulin dimer in the core of the microtubule cap (Fig. 10), we see that the distribution of GTP tubulin dimer decreases exponentially with the distance (defined Figure 7. Association rates. Association rates k z;eq at the 2 end as a function of k z;eq at the + end when both ends have the same binding energies G L and G T l . doi:10.1371/journal.pone.0006378.g007  by l) from the reference point, the base of the crown. These data generated fit with the expression 0:1e l=l except near the crown base where the distribution of GTP tubulin dimer grows faster than an exponential. Therefore we can say that parameter l is a good measure of the cap length (in tubulin dimer units). We have plotted the values of l against the relative concentration of GTP tubulin dimer, d, in Fig. 11 showing that the cap size grows linearly with d.
It is interesting to note that if one had a specific marker for GTP-tubulin it might not be easy to measure directly the length of the GTP cap. Nevertheless, one could measure how the GTPtubulin density decreases as one moves away from the tip of the microtubule. Fitting that density to an exponential curve, one would be able to measure the characteristic length l and compare it to the prediction of our model shown in Fig. 11.
Fluctuations in the length of the microtubule cap at the + end is evident over time as shown in Fig. 12 for two free tubulin dimer concentrations, d~2 and d~5. The sharp decreases in length occur when a single GTP-tubulin dimer takes a long time to hydrolyse to GDP dimer. For example, between the arrows in Fig. 12 the hydrolysis time is 30 s compared to the average GTP hydrolysis time of 3 s. The decreased length depicts the new length of the cap which is from the next closest GTP-tubulin dimer to the base of the crown. Moreover we can evaluate the length of the cap and this averages approximately 60 dimers for d~2 and 300 dimers for d~5. Another way of evaluating the size of the cap is to count the number of GTP-tubulin dimers over the 13 protofilaments. These data are shown in Fig. 13 and fluctuations in the number of GTP-tubulin dimers is apparent further supporting the variations in cap length seen in Fig. 14. Performing a similar analysis and considering only the crown, fluctuations in length occur from 1 to 7 GTP tubulin dimers when d~2.
It is also interesting to note that we have repeated these simulations where the number of protofilaments is 12 to 15 [1] and we see no change in the growing rate of the microtubule compared to the 13 protofilament microtubules used in our analysis. Moreover, when we change the monomer shift from the normal 3 tubulin monomers to 0-5 monomers at the seam, which occurs between protofilament 1 and 13, we again see no change in microtubule growing rates. These data indicate that the protofilament number and the displacement at the seam do not contribute directly to microtubule dynamics.

Discussion
Here we have generated a Monte Carlo model of the assembly/ disassembly of microtubules where we considered the rates of   attachment and detachment of the tubulin dimer at the + and 2 ends as well as the hydrolysis rate of the GTP to GDP tubulin dimer within the microtubule. Our starting premise is that microtubules are in an aqueous solution of polymerised and free tubulin dimers. We assume that this is a perfect solution in a thermodynamic context. At the critical concentration of free tubulin dimer for assembly the system is at thermodynamic equilibrium. At this equilibrium we use the thermodynamic laws [12] to relate the ratio between the attachment and detachment rate of a GTP tubulin dimer to its binding energy to the microtubule Eq. 8 .
In comparing the development of our model with that of Van Buren et al 2002, we take into account the contribution of water to the entropy of the system ( Eq. 8 ). We also take the association rate k z as a variable parameter rather than setting it to the average growing rate of microtubules. Furthermore, we have not set the GTP hydrolysis rate to zero at any point.
Our model is similar to the current model of Van Buren et al [9] with the following refinements. Firstly, we take into account the contribution of water to the entropy of the system which leads to the factor N w appearing in Eq. 8 and 11 . By treating the system as a solution, rather than a gas as in [9], we find that the longitudinal binding energies, G L and G l , are systematically 4 kT smaller than in [9]. Secondly, in our hands (B.P. and K.P. independently) the calculations for G L and G l in [9] can only be valid if G l~0 . This cannot be the case otherwise the protofilaments in the microtubules would fall apart. One solution to this problem is to take k z;eq as a separate parameter that differs from the average microtubule protoflament growing rate K z ( Eq. 16 ). However, this means that this value has to be determined from the experimental data and here we have done this. Thirdly, we have determined the values of G L and G l using a non-zero hydrolysis rate throughout microtubule polymerisation which was determined as 0:3=s=dimer from the dilution experiments in [10]. In [9] the hydolysis rate was set to zero to determine G L and G l and then set to 1=s=dimer for their simulations. However, if the G L and G l values are not adjusted when the hydrolysis rate is changed then the polymerisation rate of the microtubule will be less than the experimental value. Fourthly, for equation Eq. 12 to be valid the hydolysis rate has to be much smaller than the smallest value of k z;eq considered in [9] i.e. 2=s=dimer as explained in the derivation of equation Eq. 12 here in this paper.
Using our model we established that the dimer attachment rate is larger ( ca. 2 times) than the average microtubule growing rate which indicates that the microtubule dynamics is a complex stochastic process which depends crucially on the various configurations that the microtubule cap assumes ( i.e the shape, size and dimer constitution of the cap). We have fit all the parameters in our model to the experimental data of Walker et al [10,14]. The parameters are, the GTP tubulin dimer attachment rate at both ends, the lateral binding energies of GTP and GDP tubulin dimers at both ends and the longitudinal binding energy, which we assume to be the same at the + and the 2 ends, which were all determined from the microtubule polymerisation experiments detailed in Walker et al [14]. In addition, the hydrolysis rate was determined from the microtubule polymerisation dilution experiment described in Walker et al [10].
Once the model had been fit to the experimental data we examined the structure of the microtubule cap. We noted from the length distribution of GTP tubulin dimers and the differences in the hydrolysis time of GTP that the GTP tubulin dimer cap is long and can be up to several microns in length dependent on the free dimer concentration (d) used in the simulation exercise. In the microtubule polymerisation experiments performed by Walker et al, they determined the growing rate of individual microtubules over a period of 6-8 s. They replaced the free tubulin dimer solution with buffer and measured the time it took for the microtubules to start to depolymerise, the time delay value. They found that the time-delay increased slightly with the growing rate of the microtubule and concluded that the GTP cap was only one dimer long. However, using our model to simulate their data we can reproduce their graphical data almost exactly including the scattering around the linear fit (cf Fig. 5 with Fig. 3 in [10]). This indicates that their data are compatible with a long microtubule cap and that their fluctuations (the observed scatter) is due to the stochastic nature of the polymerisation process rather than experimental errors.  Finally, we have found that neither the number of protofilaments nor the topology of the microtubule seam have any significant effect on the microtubule dynamics. This is important given the natural variation in protofilament numbers in different biological systems.
While we have compared the predictions of our model with the in vitro experimental data, it is important to point out that the model also applies to in vivo microtubles but only after the parameters of the model have been modified to reflect the differences between the buffer used in in-vitro experiments and the cytoplasm of a cell. A buffer and the cytoplasm differ by the type of solutes and their concentrations and these differences affect the thermodynamics properties of the system. A variation in the ion concentrations is likely to affect the binding energies of the tubulin dimers, though probably not by very much, and the variation of the entropy DS, if the solution is not perfect. So the values of G L and G l would have to be altered. The hydrolysis rate r h is also likely to differ and as ions are likely to interfere with the binding of the free tubulin dimers, the attachment rate k z is likely to differ also. One must also bear in mind that different cells can have different tubulin isotypes and that this can also affect the values of the parameters of the model. So before one can apply our model to the in vivo situation all the model parameters would have to be fitted to the specific cell type.