The program is written using NetBeans, a Java IDE (Integrated Development Environment). The Java Swing JFrame Object is used as the basis of the GUI. An introduction to easily making GUI components with Netbeans is available here. The MVC (Model-View-Controller) design pattern is used to construct the user interface. Execution begins in MetabolismForm.java's main function, which creates a new MetabolismForm object and shows it. The metabolism form starts off by making several new objects, the chemical kinetic engine, two graph models, dialogue boxes and graphs and tables. The constructor initialises the components in initComponents() which is automatically written when a GUI is designed with Netbeans. The JLineGraph and JScatterGraph objects are added to a JPanel object called graphPanel (More about these later).
The chemical engine constructor randomly generates the triangular matrix of 'spins' (of up to 250 in length), with values ranging from -1 to 1. These are fixed throughout a run, and are used to determine the Gibbs Free Energy of Formation of a molecule. Several constants are defined as static variables in Engine: R, Temperature, the volume of the reactor, and a threshold rate (minRelaventRate = 0.0001). If a chemical has a rate of reaction less than minRelaventRate, that reaction is deemed negligable during the timescale of the simulation, and so is set to zero. This is done to reduce the number of reactions that need to be simulated. The effect of altering this minimum rate is explored in the control experiments section.
The addMolecule() method of engine is called twice, once to initialise 10000 molecules of type A, and again to initialise 10000 molecules of type B. addMolecule() works by passing the string of the new molecule to be added to getMolecule(). getMolecule() first checks whether this molecule already exists in the HashMap MoleculesByString. If it does not exist then a new molecule object is constructed and a reference to the new molecule object returned. When a new molecule is constructed the energy of that molecule is calculated in the method getEnergy(String), and the new molecule object is added to the moleculesByString ArrayList. If a molecule of this (formula/sequence/type) already exists then there is no need to construct a new molecule object, and a reference to the old molecule object of this type is returned instead. The method make(int) of the Molecule class is then applied to the returned Molecule object, in addMolecule().
The make(int number) method of the Molecule class first checks whether this molecule is allowed to increase in concentration (i.e. FixedCount == 0). Some molecule types can have their concentrations fixed in the flow reactor, being replenished magically at exactly the same rate at which they are used up. If FixedCount == 1, the method returns having not increased the number of molecules of this type. [N.B. It should be clear that a molecule object represents a type of molecule, not a single molecule of that type.] If FixedCount == 0, then the count of molecules of this type is increased by 'number'. If make is called on a Molecule object, which represents a molecule type which has not previously existed, i.e. if the count variable for this type has been zero previously, then a molecule instance variable 'exists' is set to true, and this molecule is added to the existingMolecules ArrayList of the Engine. Since this molecule type now has a non-zero count, it will be necessary to consider all the reactions for which this molecule type is a reactant. These are all the possible ligation reactions and cleavage reactions which this molecule type can undergo, and are calculated in methods calculateLigationReactions() and calculateCleavageReactions() respectively.
The Molecule method calculateLigationReactions() first defines a reference to a Reaction object, called reaction. Then it creates an iterator to the ArrayList existingMolecules. It uses this iterator to cycle once through all the molecules in the existingMolecules ArrayList. One cycle is as follows. The next molecule in existingMolecules is chosen and called m. The ligation (joining) reaction between 'this' molecule, and the m molecule is calculated in the Engine method getReaction(String, String), and the output of this method is a reference to a reaction, which is assigned to 'reaction' above. getReaction(String, String) takes two strings (molecular formule) and concatenates the first to the second (only in that order). It then checks to see if the HashMap 'reactions' contains this ligation reaction. The HasKey is just the sequence of the ligation product. If this reaction is contained in HashMap, then the reaction is returned to calculateLigationReactions(), otherwise PossibleReaction(String,String) is called.
PossibleReaction(String, String)'s job is to make, return and manage a new reaction object. getMolecule is called for the two reactants, and the ligation product string. Remember, getMolecule creates a new Molecule object if that molecule object was not in MoleculesByString. Remember also that when a new Molecule object is constructed, getEnergy is called. getEnergy(String) returns the Gibbs Free Energy of the molecule. This is calculated by assuming that A atoms = 1 and B atoms = -1. The spin-vector corresponding to the length of the molecule string is chosen, and the dot product of these two vectors is defined as the energy of the molecule. PossibleReaction then calculates the difference in Gibbs Free Energies (dG) between the products and the reactants, and applies the equation K = Exp(-dG/RT) to calculate the equailibrium ratio from dG, and the constant R, at temperature T. In the current implementation the forward rate is set to K, and the backward rate is set to forward_rate/ K. This results in a ratio of forward and backward rates equal to K. A reference to a reaction is defined. If either the forward rate or the backward rate is greater than minRelaventRate and if the length of the product sequence is less than an integer number e.g. 6 or 10, then a new reaction instance is created, and that reaction is assigned the forward and backward rates. This new reaction is put into the reactions HashMap straight away. However, if the above conditions are not met then a null reaction is returned by PossibleReaction to getReaction and then to calculateLigationReaction() again. The net result is that calculateLigationReaction has obtained a reaction between 'this' molecule and m. If the reaction is not null, then a new ReactionEvent is created called event, and this is added to the reactionEvents prority queue which will be used to implement the next reaction method. The molecule m keeps a record of the reactionEvents in which it takes part. In the general case where this != m, the ligation reaction is calculated for m + this, as well as for this + m.
The Molecule method calculateCleavageReactions() determines all possible pairs of substrings of molecule 'this' that can be produced by cleavage at one site. Just as in calculateLigationReactions(), getReaction is used to obtain the reaction object for the cleavage reaction, and if this is not null, an event is added to to the priority queue. Thus ends the make() function which was called from addMolecule, which was called from the constructor of MoleculeForm.
The above initialisation steps have been necessary to define the chemical network topology for the possible (moleculesByString) and actual (existingMolecules) molecules, i.e. to construct the appropriate set of reactions that need to be considered and the corresponding ReactionEvents for the priority queue. Next, the priority queue must be kept properly updated when molecule numbers change. This is done in recalculateOtherEvents.
Now, for the two types of molecule that exist, "A" and "B", recalculateOtherEvents is called, with a null ReactionEvent. This goes through all the reactionEvents stored in this Molecule and applies recalculate to each event. This must be done because the number of molecules of type A and B have changed and so the propensities of their reactions has changed. recalculate() sets the new time for all events that this molecule is involved in, based on the current propensity of the reaction, and calls reschedule to update the priority queue. recalculateOtherEvents can also be called with a non-null ReactionEvent, in which case, the propensity is not recalculated for that event. This will be necessary in the main simulation loop.
Next the JTable is passed the table model mtm. See documentation on JTable here. The cgm and mgm models are updated, and the GUI is the updated. Scroll down to jButton1ActionPerformed. This is the method called when the 'advance' button is pressed. The button works as a toggle for the Boolean variable runEngine. The user enters the number of iterations per visualization in the jSpinner1 object, and this number is stored in numEvents. Some final methods are defined. If runEngine == true, then the button text is changed to say stop, and a new SwingWorker object is created. Some SwingWorker methods are written, and the worker is started. Here an infinate while loop calls the doEvents method of engine, updates the cgm (records data) and calls updateAComponent method defined above. The actual work of the model is done in doEvents which is called with argument numEvents. This short method calls doNextEvent() numEvent times. doNextEvent gets the event that is highest in the priority queue and stores it as a ReactionEvent object called event. The products and reactent molecule objects of this ReactionEvent are obtained and stored in comb, c1 and c2. If the event is a cleavage event then the make function we looked at earlier is called for the two products and destroy is called for the reactant. Otherwise the two products are destroyed and the reactant is created. Arguments to make and destroy are 1, since only 1 molecule of each type is made and destroyed. The currentTime is set to the event time and updateEventQueue updates the priority queue. Make and destroy deals with all the complexity of maintianing a valid reaction network.