The above equations suffice for calculating the probability that a given template site will be occupied at equilibrium. However, we do not wish to give our DNA systems time to reach equilibrium, due to the competition problem discussed above. Therefore, our simulation also needs to concern itself with the dynamical behavior of DNA interactions through time--i.e.,how long will it take for an oligo to bind to a given site, and how long will it be before the oligo disassociates again?
Such dynamic information about chemical reactions is usually expressed in terms of a rate constant, which tells you what fraction of the reactants will be turned into products per unit time. This is broken down into a forward rate constant kf and a reverse rate constant kr. In our case, we will consider the forward rate constant to be the fraction of unoccupied sites that are converted to occupied sites, per unit time. To derive this, we used an empirical formula derived from Wetmur and Davidson's work [Wetmur-Davidson-68]:
Where c is again the concentration of oligos, and L is the length of the duplexes in base pairs. Note that this formula is not dependent on the equilibrium constant. It really expresses the time required for a template site to bump into an oligo randomly wandering around in solution. The reason this works is that it is thought that the rate-limiting step in the association of DNA strands is the time to make this initial contact, after which the two strand quickly "zip up" to form the complete duplex. Note that the above formula is proportional to the oligo concentration, which makes sense because the more oligos there are, the less time should be required to bump into one at random. The factor 1/sqrt(L) corresponds to the relationship between the speed of Brownian motion of an object (in this case, an oligo strand) and its mass.
The reverse rate constant kr for a reaction can be derived from kf if the equilibrium ratio of products to reactants is known. If c1 is the equilibrium concentration of reactants, and c2 is the equilibruim concentration of products, then (see this figure) we know that
so
In our case, we already know how to calculate c2/c1; it is just the ratio r of bound to unbound sites that we showed how to calculate in the previous subsection.
Figure 5 : Illustration of interconversion between reactants and products. At equilibrium the total flow to the right (from reactants to products) must equal the flow to the left (products to reactants)
This formula is used in the simulator to derive the probability that a given unnocupied site will become occupied, and the probability that an occupied site will become unoccupied (just switch kf and kr), during each tick of the simulation.
This essentially sums up the physical chemistry used in the simulator.