Skip to content
Browse files

Initial commit

  • Loading branch information...
0 parents commit 98f391ac71ffac18b226e6d8e7fcfb11a134832d @satya-arjunan satya-arjunan committed Apr 20, 2010
33 CoordinateLogProcess.cpp
@@ -0,0 +1,33 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "CoordinateLogProcess.hpp"
+
+LIBECS_DM_INIT(CoordinateLogProcess, Process);
110 CoordinateLogProcess.hpp
@@ -0,0 +1,110 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __CoordinateLogProcess_hpp
+#define __CoordinateLogProcess_hpp
+
+#include <fstream> //provides ofstream
+#include <MethodProxy.hpp>
+#include "IteratingLogProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_CLASS(CoordinateLogProcess, IteratingLogProcess)
+{
+public:
+ LIBECS_DM_OBJECT(CoordinateLogProcess, Process)
+ {
+ INHERIT_PROPERTIES(IteratingLogProcess);
+ }
+ CoordinateLogProcess():
+ theMoleculeSize(0) {}
+ virtual ~CoordinateLogProcess() {}
+ virtual void initializeLastOnce()
+ {
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ theMoleculeSize += theProcessSpecies[i]->size();
+ }
+ theLogFile.open(FileName.c_str(), ios::trunc);
+ initializeLog();
+ logSpecies();
+ }
+ virtual void fire()
+ {
+ if(theTime <= LogDuration)
+ {
+ logSpecies();
+ theTime += theStepInterval;
+ }
+ else
+ {
+ theTime = libecs::INF;
+ theLogFile.flush();
+ theLogFile.close();
+ }
+ thePriorityQueue->moveTop();
+ }
+ void logSpecies()
+ {
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ logMolecules(i);
+ }
+ theLogFile << endl;
+ }
+protected:
+ void initializeLog()
+ {
+ Point aCenterPoint(theSpatiocyteStepper->getCenterPoint());
+ theLogFile
+ << "startCoord:" << theSpatiocyteStepper->getStartCoord()
+ << " rowSize:" << theSpatiocyteStepper->getRowSize()
+ << " layerSize:" << theSpatiocyteStepper->getLayerSize()
+ << " colSize:" << theSpatiocyteStepper->getColSize()
+ << " width:" << aCenterPoint.z*2
+ << " height:" << aCenterPoint.y*2
+ << " length:" << aCenterPoint.x*2
+ << " voxelRadius:" << theSpatiocyteStepper->getVoxelRadius()
+ << " moleculeSize:" << theMoleculeSize << endl;
+ }
+ void logMolecules(int anIndex)
+ {
+ Species* aSpecies(theProcessSpecies[anIndex]);
+ for(unsigned int i(0); i != aSpecies->size(); ++i)
+ {
+ theLogFile << ", " << aSpecies->getCoord(i);
+ }
+ }
+private:
+ double theMoleculeSize;
+};
+
+#endif /* __CoordinateLogProcess_hpp */
300 DiffusionInfluencedReactionProcess.cpp
@@ -0,0 +1,300 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "DiffusionInfluencedReactionProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_INIT(DiffusionInfluencedReactionProcess, Process);
+
+
+void DiffusionInfluencedReactionProcess::initializeSecond()
+{
+ ReactionProcess::initializeSecond();
+ A->setDiffusionInfluencedReactantPair(B);
+ B->setDiffusionInfluencedReactantPair(A);
+ r_v = theSpatiocyteStepper->getVoxelRadius();
+ D_A = A->getDiffusionCoefficient();
+ D_B = B->getDiffusionCoefficient();
+ calculateReactionProbability();
+ if(A->getIsDiffusing())
+ {
+ A->setDiffusionInfluencedReaction(this, B->getID(), p);
+ }
+ if(B->getIsDiffusing())
+ {
+ B->setDiffusionInfluencedReaction(this, A->getID(), p);
+ }
+}
+
+void DiffusionInfluencedReactionProcess::initializeThird()
+{
+}
+
+//Do the reaction A + B -> C + D. So that A <- C and B <- D.
+//We need to consider that the source molecule can be either A or B.
+//If A and C belong to the same compartment, A <- C.
+//Otherwise, find a vacant adjoining voxel of A, X which is the same compartment
+//as C and X <- C.
+//Similarly, if B and D belong to the same compartment, B <- D.
+//Otherwise, find a vacant adjoining voxel of C, Y which is the same compartment
+//as D and Y <- D.
+bool DiffusionInfluencedReactionProcess::react(Voxel* moleculeB, Voxel** target)
+{
+ Voxel* moleculeA(*target);
+ //moleculeA is the source molecule. It will be soft-removed (id kept intact)
+ //by the calling Species if this method returns true.
+ //moleculeB is the target molecule, it will also be soft-removed by the
+ //calling Species.
+ //First let us make sure moleculeA and moleculeB belong to the
+ //correct species.
+ if(moleculeA->id != A->getID())
+ {
+ Voxel* tempA(moleculeA);
+ moleculeA = moleculeB;
+ moleculeB = tempA;
+ }
+ //If the product C is not in the same compartment as A,
+ //we need to find a vacant adjoining voxel of A that belongs
+ //to the compartment of C:
+ Voxel* moleculeC;
+ if(A->getVacantID() != C->getVacantID())
+ {
+ moleculeC = C->getRandomAdjoiningVoxel(moleculeA);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of A which can be occupied by C:
+ if(moleculeC == NULL)
+ {
+ return false;
+ }
+ }
+ else
+ {
+ moleculeC = moleculeA;
+ }
+ //If it has two products:
+ if(D != NULL)
+ {
+ //If the product D is not in the same compartment as B,
+ //we need to find a vacant adjoining voxel of C that belongs
+ //to the compartment of D:
+ Voxel* moleculeD;
+ if(B->getVacantID() != D->getVacantID())
+ {
+ moleculeD = D->getRandomAdjoiningVoxel(moleculeC, moleculeC);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of A which can be occupied by C:
+ if(moleculeD == NULL)
+ {
+ return false;
+ }
+ }
+ else
+ {
+ moleculeD = moleculeB;
+ }
+ if(E != NULL)
+ {
+ //moleculeA and moleculeB will not be selected here because
+ //they are still not vacant:
+ Voxel* moleculeE(E->getRandomAdjoiningVoxel(moleculeD, moleculeC,
+ moleculeD));
+ if(moleculeE == NULL)
+ {
+ return false;
+ }
+ E->addMolecule(moleculeE);
+ }
+ else if(variableE != NULL)
+ {
+ variableE->addValue(1);
+ }
+ moleculeB->id = B->getVacantID();
+ D->addMolecule(moleculeD);
+ }
+ else
+ {
+ //Hard remove the B molecule since this is a single product reaction:
+ moleculeB->id = B->getVacantID();
+ }
+ //Hard remove the A molecule, in case C is in a different compartment:
+ moleculeA->id = A->getVacantID();
+ C->addMolecule(moleculeC);
+ return true;
+}
+
+void DiffusionInfluencedReactionProcess::finalizeReaction()
+{
+ //The number of molecules may have changed for both reactant and product
+ //species. We need to update SpatiocyteNextReactionProcesses which are
+ //dependent on these species:
+ for(vector<ReactionProcess*>::const_iterator
+ i(theInterruptingProcesses.begin());
+ i!=theInterruptingProcesses.end(); ++i)
+ {
+ (*i)->substrateValueChanged(theSpatiocyteStepper->getCurrentTime());
+ }
+}
+
+void DiffusionInfluencedReactionProcess::calculateReactionProbability()
+{
+ //Refer to the paper for the description of the variables used in this
+ //method.
+ if(A->getIsVolume() && B->getIsVolume())
+ {
+ if(A != B)
+ {
+ if(p == -1)
+ {
+ p = k/(6*sqrt(2)*(D_A+D_B)*r_v);
+ }
+ else
+ {
+ k = p*(6*sqrt(2)*(D_A+D_B)*r_v);
+ }
+ }
+ else
+ {
+ if(p == -1)
+ {
+ p = k/(6*sqrt(2)*D_A*r_v);
+ }
+ else
+ {
+ k = p*(6*sqrt(2)*D_A*r_v);
+ }
+ }
+ }
+ else if(!A->getIsVolume() && !B->getIsVolume())
+ {
+ if(A != B)
+ {
+ if(p == -1)
+ {
+ p = pow(2*sqrt(2)+4*sqrt(3)+3*sqrt(6)+sqrt(22), 2)*k/
+ (72*(6*sqrt(2)+4*sqrt(3)+3*sqrt(6))*(D_A+D_B));
+ }
+ else
+ {
+ k = p*(72*(6*sqrt(2)+4*sqrt(3)+3*sqrt(6))*(D_A+D_B))/
+ pow(2*sqrt(2)+4*sqrt(3)+3*sqrt(6)+sqrt(22), 2);
+ }
+ }
+ else
+ {
+ if(p == -1)
+ {
+ p = pow(2*sqrt(2)+4*sqrt(3)+3*sqrt(6)+sqrt(22), 2)*k/
+ (72*(6*sqrt(2)+4*sqrt(3)+3*sqrt(6))*(D_A));
+ }
+ else
+ {
+ k = p*(72*(6*sqrt(2)+4*sqrt(3)+3*sqrt(6))*(D_A))/
+ pow(2*sqrt(2)+4*sqrt(3)+3*sqrt(6)+sqrt(22), 2);
+ }
+ }
+ }
+ else if(A->getIsVolume() && B->getIsLipid())
+ {
+ if(p == -1)
+ {
+ p = 24*k*r_v/((6+3*sqrt(3)+2*sqrt(6))*D_A);
+ }
+ else
+ {
+ k = p*((6+3*sqrt(3)+2*sqrt(6))*D_A)/(24*r_v);
+ }
+ }
+ else if(A->getIsLipid() && B->getIsVolume())
+ {
+ if(p == -1)
+ {
+ p = 24*k*r_v/((6+3*sqrt(3)+2*sqrt(6))*D_B);
+ }
+ else
+ {
+ k = p*((6+3*sqrt(3)+2*sqrt(6))*D_B)/(24*r_v);
+ }
+ }
+ else if(A->getIsVolume() && !B->getIsVolume())
+ {
+ if(p == -1)
+ {
+ p = sqrt(2)*k/(3*D_A*r_v);
+ }
+ else
+ {
+ k = p*(3*D_A*r_v)/sqrt(2);
+ }
+ }
+ else if(!A->getIsVolume() && B->getIsVolume())
+ {
+ if(p == -1)
+ {
+ p = sqrt(2)*k/(3*D_B*r_v);
+ }
+ else
+ {
+ k = p*(3*D_B*r_v)/sqrt(2);
+ }
+ }
+ else
+ {
+ THROW_EXCEPTION(ValueError,
+ String(getPropertyInterface().getClassName()) +
+ "[" + getFullID().asString() +
+ "]: Error in type of second order reaction.");
+ }
+}
+
+void DiffusionInfluencedReactionProcess::printParameters()
+{
+ cout << "[" <<
+ A->getVariable()->getSystemPath().asString() << ":" <<
+ A->getVariable()->getID() << "] + [" <<
+ B->getVariable()->getSystemPath().asString() << ":" <<
+ B->getVariable()->getID() << "] -> [" <<
+ C->getVariable()->getSystemPath().asString() << ":" <<
+ C->getVariable()->getID() << "]";
+ if(D != NULL)
+ {
+ cout << " + [" <<
+ D->getVariable()->getSystemPath().asString() << ":" <<
+ D->getVariable()->getID() << "]";
+ }
+ cout << ": k=" << k << ", p=" << p <<
+ ", p_A=" << A->getReactionProbability(B->getID()) <<
+ ", p_B=" << B->getReactionProbability(A->getID()) << endl;
+}
+
+
+
+
+
+
74 DiffusionInfluencedReactionProcess.hpp
@@ -0,0 +1,74 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __DiffusionInfluencedReactionProcess_hpp
+#define __DiffusionInfluencedReactionProcess_hpp
+
+#include "ReactionProcess.hpp"
+
+LIBECS_DM_CLASS(DiffusionInfluencedReactionProcess, ReactionProcess)
+{
+public:
+ LIBECS_DM_OBJECT(DiffusionInfluencedReactionProcess, Process)
+ {
+ INHERIT_PROPERTIES(ReactionProcess);
+ }
+ DiffusionInfluencedReactionProcess() {}
+ virtual ~DiffusionInfluencedReactionProcess() {}
+ virtual void addSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule) {}
+ virtual void removeSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule) {}
+ virtual void initialize()
+ {
+ ReactionProcess::initialize();
+ if(getOrder() != 2)
+ {
+ THROW_EXCEPTION(ValueError,
+ String(getPropertyInterface().getClassName()) +
+ "[" + getFullID().asString() +
+ "]: Only second order scheme is allowed for "+
+ "diffusion influenced reactions.");
+ }
+ }
+ virtual void initializeSecond();
+ virtual void initializeThird();
+ virtual bool react(Voxel*, Voxel**);
+ virtual void printParameters();
+ virtual void finalizeReaction();
+protected:
+ void calculateReactionProbability();
+protected:
+ double D_A;
+ double D_B;
+ double r_v;
+ double V;
+};
+
+#endif /* __DiffusionInfluencedReactionProcess_hpp */
33 DiffusionProcess.cpp
@@ -0,0 +1,33 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "DiffusionProcess.hpp"
+
+LIBECS_DM_INIT(DiffusionProcess, Process);
256 DiffusionProcess.hpp
@@ -0,0 +1,256 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __DiffusionProcess_hpp
+#define __DiffusionProcess_hpp
+
+#include <sstream>
+#include <MethodProxy.hpp>
+#include "SpatiocyteProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_CLASS(DiffusionProcess, SpatiocyteProcess)
+{
+ typedef const void (DiffusionProcess::*WalkMethod)(void) const;
+public:
+ LIBECS_DM_OBJECT(DiffusionProcess, Process)
+ {
+ INHERIT_PROPERTIES(Process);
+ PROPERTYSLOT_SET_GET(Real, D);
+ PROPERTYSLOT_SET_GET(Real, P);
+ PROPERTYSLOT_SET_GET(Real, WalkProbability);
+ }
+ DiffusionProcess():
+ D(0),
+ P(1),
+ WalkProbability(1),
+ theWalkMethod(&DiffusionProcess::volumeWalk) {}
+ virtual ~DiffusionProcess() {}
+ SIMPLE_SET_GET_METHOD(Real, D);
+ SIMPLE_SET_GET_METHOD(Real, P);
+ SIMPLE_SET_GET_METHOD(Real, WalkProbability);
+ virtual void initialize()
+ {
+ SpatiocyteProcess::initialize();
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->setDiffusionCoefficient(D);
+ }
+ }
+ virtual void initializeThird()
+ {
+ Species* aSpecies(theProcessSpecies[0]);
+ isVolume = aSpecies->getIsVolume();
+ double rho(aSpecies->getMaxReactionProbability());
+ if(D > 0)
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ if((*i)->getIsVolume() != isVolume)
+ {
+ THROW_EXCEPTION(ValueError, String(
+ getPropertyInterface().getClassName()) +
+ "[" + getFullID().asString() +
+ "]: A DiffusionProcess can only execute" +
+ " multiple species when they are all either" +
+ " in a volume compartment or a surface" +
+ " compartment, not both concurrently. " +
+ getIDString(theProcessSpecies[0]) + " and " +
+ getIDString(*i) + " belong to different" +
+ " types of compartment.");
+ }
+ if(rho < (*i)->getMaxReactionProbability())
+ {
+ if(rho > P)
+ {
+ THROW_EXCEPTION(ValueError, String(
+ getPropertyInterface().getClassName()) +
+ "[" + getFullID().asString() +
+ "]: Create separate" +
+ " DiffusionProcesses for " +
+ getIDString(aSpecies) + " and " +
+ getIDString(*i) + " since their" +
+ " reaction probabilities are not the" +
+ " same and the latter's reaction" +
+ " probability is higher than P.");
+ }
+ aSpecies = *i;
+ rho = (*i)->getMaxReactionProbability();
+ }
+ }
+ }
+ if(rho > P)
+ {
+ WalkProbability = P/rho;
+ }
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->rescaleReactionProbabilities(WalkProbability);
+ }
+ if(D > 0)
+ {
+ double r_v(theSpatiocyteStepper->getVoxelRadius());
+ double lambda(2.0/3);
+ if(!isVolume)
+ {
+ lambda = pow((2*sqrt(2)+4*sqrt(3)+3*sqrt(6)+sqrt(22))/
+ (6*sqrt(2)+4*sqrt(3)+3*sqrt(6)), 2);
+ }
+ theStepInterval = lambda*r_v*r_v*WalkProbability/D;
+ }
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->setDiffusionInterval(theStepInterval);
+ }
+ if(isVolume)
+ {
+ if(WalkProbability == 1)
+ {
+ theWalkMethod = &DiffusionProcess::volumeWalk;
+ }
+ else
+ {
+ theWalkMethod = &DiffusionProcess::volumeWalkCollide;
+ }
+ }
+ else
+ {
+ if(WalkProbability == 1)
+ {
+ theWalkMethod = &DiffusionProcess::surfaceWalk;
+ }
+ else
+ {
+ theWalkMethod = &DiffusionProcess::surfaceWalkCollide;
+ }
+ }
+ //At the start of the simulation, we must make sure the CollisionProcess
+ //is fired first before the DiffusionProcess. This is to make sure
+ //the reaction probability is valid for reactants that are initially
+ //at contact:
+ }
+ virtual void printParameters()
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ cout << getIDString(*i) << " ";
+ }
+ cout << ":" << endl << " Diffusion interval=" << theStepInterval <<
+ ", D=" << D << ", Walk probability (P/rho)=" <<
+ WalkProbability << endl;
+ }
+ virtual void fire()
+ {
+ (this->*theWalkMethod)();
+ theTime += theStepInterval;
+ thePriorityQueue->moveTop();
+ }
+ const void volumeWalk() const
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->volumeWalk();
+ }
+ }
+ const void volumeWalkCollide() const
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->volumeWalkCollide();
+ }
+ }
+ const void surfaceWalk() const
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->surfaceWalk();
+ }
+ }
+ const void surfaceWalkCollide() const
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->surfaceWalkCollide();
+ }
+ }
+ virtual void initializeLastOnce()
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->addInterruptedProcess(
+ reinterpret_cast<SpatiocyteProcess*>(this));
+ }
+ }
+ virtual void addSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule)
+ {
+ if(theStepInterval == libecs::INF)
+ {
+ theStepInterval = theProcessSpecies[0]->getDiffusionInterval();
+ theTime = theSpatiocyteStepper->getCurrentTime() + theStepInterval;
+ thePriorityQueue->move(theQueueID);
+ }
+ }
+ virtual void removeSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule)
+ {
+ if(theStepInterval != libecs::INF)
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ if((*i)->size())
+ {
+ return;
+ }
+ }
+ theStepInterval = libecs::INF;
+ theTime = theStepInterval;
+ thePriorityQueue->move(theQueueID);
+ }
+ }
+protected:
+ bool isVolume;
+ double D;
+ double P;
+ double WalkProbability;
+ WalkMethod theWalkMethod;
+};
+
+#endif /* __DiffusionProcess_hpp */
102 FluorescentProteinImagingProcess.cpp
@@ -0,0 +1,102 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "FluorescentProteinImagingProcess.hpp"
+
+LIBECS_DM_INIT(FluorescentProteinImagingProcess, Process);
+
+void FluorescentProteinImagingProcess::incSpeciesLatticeCount()
+{
+ for(unsigned int i(0); i != thePositiveSpecies.size(); ++i)
+ {
+ Species* aSpecies(thePositiveSpecies[i]);
+ unsigned int aMoleculeSize(aSpecies->size());
+ for(unsigned int j(0); j != aMoleculeSize; ++j)
+ {
+ Voxel* aMolecule(aSpecies->getMolecule(j));
+ ++theLattice[theProcessSpeciesIndices[i][0]][
+ aMolecule->coord-theStartCoord];
+ for(unsigned int k(1); k != theProcessSpeciesIndices[i].size(); ++k)
+ {
+ Voxel* anAdjoin(aMolecule->adjoiningVoxels[k-1]);
+ ++theLattice[theProcessSpeciesIndices[i][k]][
+ anAdjoin->coord-theStartCoord];
+ }
+ }
+ }
+}
+
+void FluorescentProteinImagingProcess::logFluorescentSpecies()
+{
+ vector<int> coordList;
+ for(unsigned int i(0); i != theLatticeSize; ++i)
+ {
+ for(unsigned int j(0); j != theProcessSpecies.size(); ++j)
+ {
+ if(theLattice[j][i])
+ {
+ coordList.push_back(i);
+ break;
+ }
+ }
+ }
+ int aDataSize(0);
+ streampos aStartPos(theLogFile.tellp());
+ // write the next size (create a temporary space for it)
+ theLogFile.write((char*)(&aDataSize), sizeof(aDataSize));
+ double aCurrentTime(theSpatiocyteStepper->getCurrentTime());
+ theLogFile.write((char*)(&aCurrentTime), sizeof(aCurrentTime));
+ unsigned int coordSize(coordList.size());
+ theLogFile.write((char*)(&coordSize), sizeof(coordSize));
+ for(unsigned int i(0); i != coordSize; ++i)
+ {
+ unsigned int aCoord(coordList[i]+theStartCoord);
+ theLogFile.write((char*)(&aCoord), sizeof(aCoord));
+ }
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ for(unsigned int j(0); j != coordSize; ++j)
+ {
+ unsigned int frequency(theLattice[i][coordList[j]]);
+ theLogFile.write((char*)(&frequency), sizeof(frequency));
+ }
+ }
+ aDataSize = (theLogFile.tellp()-aStartPos)-sizeof(aDataSize);
+ int aPrevDataSize(theLogFile.tellp()-theStepStartPos+sizeof(int)*2);
+ theStepStartPos = theLogFile.tellp();
+ // write the prev size at the end of this step
+ theLogFile.write((char*)(&aPrevDataSize), sizeof(aPrevDataSize));
+ streampos aCurrentPos(theLogFile.tellp());
+ theLogFile.seekp(aStartPos);
+ // write the next size at the beginning of this step
+ theLogFile.write((char*) (&aDataSize), sizeof(aDataSize));
+ theLogFile.seekp(aCurrentPos);
+}
+
185 FluorescentProteinImagingProcess.hpp
@@ -0,0 +1,185 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __FluorescentProteinImagingProcess_hpp
+#define __FluorescentProteinImagingProcess_hpp
+
+#include "VisualizationLogProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_CLASS(FluorescentProteinImagingProcess, VisualizationLogProcess)
+{
+public:
+ LIBECS_DM_OBJECT(FluorescentProteinImagingProcess, Process)
+ {
+ INHERIT_PROPERTIES(VisualizationLogProcess);
+ PROPERTYSLOT_SET_GET(Integer, MeanCount);
+ PROPERTYSLOT_SET_GET(Real, ExposureTime);
+ }
+ FluorescentProteinImagingProcess():
+ MeanCount(-1),
+ theLastExposedTime(0) {}
+ virtual ~FluorescentProteinImagingProcess() {}
+ SIMPLE_SET_GET_METHOD(Integer, MeanCount);
+ SIMPLE_SET_GET_METHOD(Real, ExposureTime);
+ virtual void initializeFourth()
+ {
+ theLatticeSize = theSpatiocyteStepper->getLatticeSize();
+ theStartCoord = theSpatiocyteStepper->getStartCoord();
+ theProcessSpecies.resize(0);
+ for(VariableReferenceVector::iterator
+ i(theNegativeVariableReferences.begin());
+ i != theNegativeVariableReferences.end(); ++i)
+ {
+ Variable* aVariable((*i).getVariable());
+ Species* aSpecies(theSpatiocyteStepper->getSpecies(aVariable));
+ if(!aSpecies->getIsLipid() &&
+ find(theProcessSpecies.begin(), theProcessSpecies.end(),
+ aSpecies) == theProcessSpecies.end())
+ {
+ theProcessSpecies.push_back(aSpecies);
+ }
+ }
+ VariableReferenceVector::iterator aNegativeVariableIter(
+ theNegativeVariableReferences.begin());
+ for(VariableReferenceVector::iterator
+ i(thePositiveVariableReferences.begin());
+ i != thePositiveVariableReferences.end(); ++i)
+ {
+ int aPositiveCoefficient((*i).getCoefficient());
+ Variable* aVariable((*i).getVariable());
+ Species* aSpecies(theSpatiocyteStepper->getSpecies(aVariable));
+ if(!aSpecies->getIsLipid())
+ {
+ if(aPositiveCoefficient > 0)
+ {
+ thePositiveSpecies.push_back(aSpecies);
+ vector<int> aProcessSpeciesIndices;
+ do
+ {
+ int aNegativeCoefficient(
+ (*aNegativeVariableIter).getCoefficient());
+ Variable* aNegativeVariable(
+ (*aNegativeVariableIter).getVariable());
+ Species* aNegativeSpecies(
+ theSpatiocyteStepper->getSpecies(aNegativeVariable));
+ int aProcessSpeciesIndex(
+ find(theProcessSpecies.begin(), theProcessSpecies.end(),
+ aNegativeSpecies) - theProcessSpecies.begin());
+ while(aNegativeCoefficient)
+ {
+ aNegativeCoefficient += 1;
+ aPositiveCoefficient -= 1;
+ aProcessSpeciesIndices.push_back(
+ aProcessSpeciesIndex);
+ }
+ ++aNegativeVariableIter;
+ }
+ while(aPositiveCoefficient);
+ theProcessSpeciesIndices.push_back(aProcessSpeciesIndices);
+ }
+ }
+ }
+ theLattice.resize(theProcessSpecies.size());
+ for(unsigned int i(0); i != theLattice.size(); ++i)
+ {
+ theLattice[i].resize(theLatticeSize);
+ }
+ resetLattice();
+ if(MeanCount != -1)
+ {
+ theStepInterval = ExposureTime/MeanCount;
+ }
+ else
+ {
+ for(vector<Species*>::const_iterator i(thePositiveSpecies.begin());
+ i != thePositiveSpecies.end(); ++i)
+ {
+ if((*i)->getDiffusionInterval() < theStepInterval)
+ {
+ theStepInterval = (*i)->getDiffusionInterval();
+ }
+ Species* reactantPair((*i)->getDiffusionInfluencedReactantPair());
+ if(reactantPair != NULL &&
+ reactantPair->getDiffusionInterval() < theStepInterval)
+ {
+ theStepInterval = reactantPair->getDiffusionInterval();
+ }
+ }
+ MeanCount = (int)ceil(ExposureTime/theStepInterval);
+ }
+ theMeanCount = (unsigned int)MeanCount;
+ }
+ virtual void initializeLastOnce()
+ {
+ ostringstream aFilename;
+ aFilename << FileName << ends;
+ theLogFile.open(aFilename.str().c_str(), ios::binary | ios::trunc);
+ initializeLog();
+ logSurfaceVoxels();
+ }
+ virtual void fire()
+ {
+ incSpeciesLatticeCount();
+ if(theTime-theLastExposedTime >= ExposureTime)
+ {
+ theLastExposedTime = theTime;
+ logFluorescentSpecies();
+ resetLattice();
+ }
+ theTime += theStepInterval;
+ thePriorityQueue->moveTop();
+ }
+protected:
+ void incSpeciesLatticeCount();
+ void logFluorescentSpecies();
+ void resetLattice()
+ {
+ for(unsigned int i(0); i != theLattice.size(); ++i)
+ {
+ for(unsigned int j(0); j != theLatticeSize; ++j)
+ {
+ theLattice[i][j] = 0;
+ }
+ }
+ }
+protected:
+ unsigned int theStartCoord;
+ unsigned int theLatticeSize;
+ int MeanCount;
+ double ExposureTime;
+ double theLastExposedTime;
+ vector<Species*> thePositiveSpecies;
+ vector<vector<int> > theLattice;
+ vector<vector<int> > theProcessSpeciesIndices;
+};
+
+#endif /* __FluorescentProteinImagingProcess_hpp */
33 IteratingLogProcess.cpp
@@ -0,0 +1,33 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "IteratingLogProcess.hpp"
+
+LIBECS_DM_INIT(IteratingLogProcess, Process);
304 IteratingLogProcess.hpp
@@ -0,0 +1,304 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __IteratingLogProcess_hpp
+#define __IteratingLogProcess_hpp
+
+#include <fstream> //provides ofstream
+#include <math.h>
+#include "SpatiocyteProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_CLASS(IteratingLogProcess, SpatiocyteProcess)
+{
+public:
+ LIBECS_DM_OBJECT(IteratingLogProcess, Process)
+ {
+ INHERIT_PROPERTIES(Process);
+ PROPERTYSLOT_SET_GET(Real, LogDuration);
+ PROPERTYSLOT_SET_GET(Real, LogInterval);
+ PROPERTYSLOT_SET_GET(Integer, Iterations);
+ PROPERTYSLOT_SET_GET(Integer, SaveInterval);
+ PROPERTYSLOT_SET_GET(Integer, Centered);
+ PROPERTYSLOT_SET_GET(Integer, InContact);
+ PROPERTYSLOT_SET_GET(Integer, Survival);
+ PROPERTYSLOT_SET_GET(Integer, RebindTime);
+ PROPERTYSLOT_SET_GET(Integer, Displacement);
+ PROPERTYSLOT_SET_GET(Integer, Diffusion);
+ PROPERTYSLOT_SET_GET(String, FileName);
+ }
+ SIMPLE_SET_GET_METHOD(Real, LogDuration);
+ SIMPLE_SET_GET_METHOD(Real, LogInterval);
+ SIMPLE_SET_GET_METHOD(Integer, Iterations);
+ SIMPLE_SET_GET_METHOD(Integer, SaveInterval);
+ SIMPLE_SET_GET_METHOD(Integer, Centered);
+ SIMPLE_SET_GET_METHOD(Integer, InContact);
+ SIMPLE_SET_GET_METHOD(Integer, Survival);
+ SIMPLE_SET_GET_METHOD(Integer, RebindTime);
+ SIMPLE_SET_GET_METHOD(Integer, Displacement);
+ SIMPLE_SET_GET_METHOD(Integer, Diffusion);
+ SIMPLE_SET_GET_METHOD(String, FileName);
+ IteratingLogProcess():
+ SpatiocyteProcess(false),
+ Centered(0),
+ Diffusion(0),
+ Displacement(0),
+ InContact(0),
+ RebindTime(0),
+ SaveInterval(-1),
+ Survival(0),
+ LogInterval(-1) {}
+ virtual ~IteratingLogProcess() {}
+ virtual void initializeSecond()
+ {
+ SpatiocyteProcess::initializeSecond();
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ if(InContact)
+ {
+ theProcessSpecies[i]->setIsInContact();
+ }
+ if(Centered)
+ {
+ theProcessSpecies[i]->setIsCentered();
+ }
+ }
+ theLogCnt = 0;
+ }
+ virtual void initializeFourth()
+ {
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ if((*i)->getDiffusionInterval() < theStepInterval)
+ {
+ theStepInterval = (*i)->getDiffusionInterval();
+ }
+ Species* reactantPair((*i)->getDiffusionInfluencedReactantPair());
+ if(reactantPair != NULL &&
+ reactantPair->getDiffusionInterval() < theStepInterval)
+ {
+ theStepInterval = reactantPair->getDiffusionInterval();
+ }
+ }
+ if(LogInterval != -1)
+ {
+ theStepInterval = LogInterval;
+ }
+ else
+ {
+ LogInterval = theStepInterval;
+ }
+ theTime = LogInterval;
+ thePriorityQueue->move(theQueueID);
+ }
+ virtual void initializeLastOnce()
+ {
+ theLogFile.open(FileName.c_str(), ios::trunc);
+ theTotalIterations = Iterations;
+ theLogValues.resize(theProcessSpecies.size());
+ if(SaveInterval == -1)
+ {
+ SaveInterval = Iterations/3;
+ }
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ if(RebindTime)
+ {
+ theLogValues[i].resize(Iterations);
+ }
+ else
+ {
+ theLogValues[i].resize((int)ceil(LogDuration/theStepInterval));
+ }
+ for(unsigned int j(0); j != theLogValues[i].size(); ++j)
+ {
+ theLogValues[i][j] = 0;
+ }
+ }
+ }
+ virtual void fire()
+ {
+ if(Iterations == 0)
+ {
+ cout << "Saving data in: " << FileName.c_str() << endl;
+ double aTime(LogInterval);
+ for(unsigned int i(0); i != theLogValues[0].size(); ++i)
+ {
+ theLogFile << setprecision(15) << aTime;
+ for(unsigned int j(0); j != theProcessSpecies.size(); ++j)
+ {
+ if(RebindTime)
+ {
+ theLogFile << "," << setprecision(15) << theLogValues[j][i];
+ }
+ else
+ {
+ theLogFile << "," << setprecision(15) << theLogValues[j][i]/
+ theTotalIterations;
+ }
+ }
+ theLogFile << endl;
+ aTime += LogInterval;
+ }
+ theLogFile.close();
+ theStepInterval = libecs::INF;
+ }
+ else if(theTime >= LogDuration)
+ {
+ theStepInterval = LogInterval;
+ theSpatiocyteStepper->reset(Iterations);
+ --Iterations;
+ cout << "Iterations left:" << Iterations << endl;
+ if(SaveInterval == 0 || Iterations%SaveInterval == 0)
+ {
+ string aFileName(FileName.c_str());
+ aFileName = aFileName + ".back";
+ cout << "Saving temporary backup data in: " << aFileName << endl;
+ ofstream aFile;
+ aFile.open(aFileName.c_str(), ios::trunc);
+ double aTime(LogInterval);
+ int completedIterations(theTotalIterations-Iterations);
+ unsigned int aSize(theLogValues[0].size());
+ if(RebindTime)
+ {
+ aSize = completedIterations;
+ }
+ for(unsigned int i(0); i != aSize; ++i)
+ {
+ aFile << setprecision(15) << aTime;
+ for(unsigned int j(0); j != theProcessSpecies.size(); ++j)
+ {
+ if(RebindTime)
+ {
+ aFile << "," << setprecision(15) <<
+ theLogValues[j][i];
+ }
+ else
+ {
+ aFile << "," << setprecision(15) <<
+ theLogValues[j][i]/completedIterations;
+ }
+ }
+ aFile << endl;
+ aTime += theStepInterval;
+ }
+ aFile.close();
+ }
+ return;
+ }
+ else
+ {
+ theSurvivalCnt = 0;
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ if(RebindTime)
+ {
+ if(theProcessSpecies[i]->getVariable()->getValue())
+ {
+ ++theSurvivalCnt;
+ }
+ if(!theLogValues[i][theTotalIterations-Iterations])
+ {
+ if(!theProcessSpecies[i]->getVariable()->getValue())
+ {
+ theLogValues[i][theTotalIterations-Iterations] =
+ theTime;
+ }
+ }
+ }
+ else if(Survival)
+ {
+ if(theProcessSpecies[i]->getVariable()->getValue())
+ {
+ ++theSurvivalCnt;
+ }
+ theLogValues[i][theLogCnt] +=
+ theProcessSpecies[i]->getVariable()->getValue()/
+ theProcessSpecies[i]->getInitMoleculeSize();
+ }
+ else if(Displacement)
+ {
+ theLogValues[i][theLogCnt] +=
+ theProcessSpecies[i]->getMeanSquaredDisplacement();
+ }
+ else if(Diffusion)
+ {
+ if(theProcessSpecies[i]->getIsVolume())
+ {
+ theLogValues[i][theLogCnt] +=
+ theProcessSpecies[i]->getMeanSquaredDisplacement()/
+ (6*theTime);
+ }
+ else
+ {
+ theLogValues[i][theLogCnt] +=
+ theProcessSpecies[i]->getMeanSquaredDisplacement()/
+ (4*theTime);
+ }
+ }
+ //By default log the values:
+ else
+ {
+ theLogValues[i][theLogCnt] +=
+ theProcessSpecies[i]->getVariable()->getValue();
+ }
+ }
+ //If all survival species are dead, go on to the next iteration:
+ if((Survival || RebindTime) && !theSurvivalCnt)
+ {
+ theStepInterval = LogDuration-theTime;
+ }
+ ++theLogCnt;
+ }
+ theTime += theStepInterval;
+ thePriorityQueue->moveTop();
+ }
+protected:
+ int Centered;
+ int Diffusion;
+ int Displacement;
+ int InContact;
+ int Iterations;
+ int RebindTime;
+ int SaveInterval;
+ int Survival;
+ int theLogCnt;
+ int theSurvivalCnt;
+ int theTotalIterations;
+ double LogDuration;
+ double LogInterval;
+ String FileName;
+ ofstream theLogFile;
+ vector<vector<double> > theLogValues;
+};
+
+#endif /* __IteratingLogProcess_hpp */
80 MassActionProcess.cpp
@@ -0,0 +1,80 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "MassActionProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_INIT(MassActionProcess, Process);
+
+void MassActionProcess::fire()
+{
+ if(theSpace == 0)
+ {
+ Species* aSpecies(NULL);
+ for(VariableReferenceVector::iterator
+ i(theVariableReferenceVector.begin());
+ i != theVariableReferenceVector.end(); ++i)
+ {
+ Variable* aVariable((*i).getVariable());
+ aSpecies = theSpatiocyteStepper->getSpecies(aVariable);
+ if(aSpecies != NULL)
+ {
+ break;
+ }
+ }
+ if(aSpecies->getIsVolume())
+ {
+ theSpace = aSpecies->getCompartment()->volume;
+ cout << "Mass Action Volume:" << theSpace << endl;
+ }
+ else
+ {
+ theSpace = aSpecies->getCompartment()->area;
+ cout << "Mass Action Area:" << theSpace << endl;
+ }
+ }
+ double velocity(k);
+ velocity *= theSpace;
+ for(VariableReferenceVector::iterator
+ s(theVariableReferenceVector.begin());
+ s != theZeroVariableReferenceIterator; ++s)
+ {
+ VariableReference aVariableReference(*s);
+ Integer aCoefficient(aVariableReference.getCoefficient());
+ do
+ {
+ ++aCoefficient;
+ velocity *= aVariableReference.getVariable()->getValue()/theSpace;
+ }
+ while(aCoefficient != 0);
+ }
+ setFlux(velocity);
+}
+
67 MassActionProcess.hpp
@@ -0,0 +1,67 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __MassActionProcess_hpp
+#define __MassActionProcess_hpp
+
+#include <ContinuousProcess.hpp>
+#include "SpatiocyteStepper.hpp"
+#include "SpatiocyteCommon.hpp"
+
+LIBECS_DM_CLASS(MassActionProcess, ContinuousProcess)
+{
+public:
+ LIBECS_DM_OBJECT(MassActionProcess, Process)
+ {
+ INHERIT_PROPERTIES(Process);
+ PROPERTYSLOT_SET_GET(Real, k);
+ }
+ MassActionProcess():
+ k(0),
+ theSpace(0) {}
+ virtual ~MassActionProcess() {}
+ SIMPLE_SET_GET_METHOD(Real, k);
+ virtual void fire();
+ virtual void initialize()
+ {
+ Process::initialize();
+ declareUnidirectional();
+ theSpatiocyteStepper =
+ reinterpret_cast<SpatiocyteStepper*>(getSuperSystem()->getStepper());
+ }
+protected:
+ double k;
+ double theSpace;
+ SpatiocyteStepper* theSpatiocyteStepper;
+};
+
+#endif /* __MassActionProcess_hpp */
+
182 OscillationAnalysisProcess.cpp
@@ -0,0 +1,182 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "OscillationAnalysisProcess.hpp"
+
+LIBECS_DM_INIT(OscillationAnalysisProcess, Process);
+
+void OscillationAnalysisProcess::testMembraneBinding()
+{
+ if(minD->size() < 100)
+ {
+ Status = 1;
+ }
+ cout << "theTime:" << theTime << " Status:" << Status << endl;
+}
+
+void OscillationAnalysisProcess::testLocalization(int aStatus)
+{
+ int quad1(0);
+ int quad2(0);
+ int quad3(0);
+ int quad4(0);
+ int aSize(minD_m->size());
+ for(int i(0); i != aSize; ++i)
+ {
+ Point aPoint(theSpatiocyteStepper->coord2point(minD_m->getCoord(i)));
+ if(aPoint.x < theCenterPoint.x/2)
+ {
+ ++quad1;
+ }
+ else if(aPoint.x >= theCenterPoint.x/2 && aPoint.x < theCenterPoint.x)
+ {
+ ++quad2;
+ }
+ else if(aPoint.x >= theCenterPoint.x && aPoint.x < theCenterPoint.x*1.5)
+ {
+ ++quad3;
+ }
+ else
+ {
+ ++quad4;
+ }
+ }
+ double threshold(0.5*aSize/4);
+ cout << "thresh:" << threshold << " size:" << aSize << endl;
+ if(quad1 < threshold || quad2 < threshold || quad3 < threshold ||
+ quad4 < threshold)
+ {
+ Status = aStatus;
+ }
+ cout << "theTime:" << theTime << " quad1:" << quad1 << " quad2:" << quad2
+ << " quad3:" << quad3 << " quad4:" << quad4 << " status:" <<
+ Status << endl;
+}
+
+void OscillationAnalysisProcess::testOscillation()
+{
+ int aSize(minD_m->size());
+ vector<double> leftPositions;
+ vector<double> rightPositions;
+ for(int i(0); i != aSize; ++i)
+ {
+ Point aPoint(theSpatiocyteStepper->coord2point(minD_m->getCoord(i)));
+ if(aPoint.x < theCenterPoint.x)
+ {
+ leftPositions.push_back(aPoint.x);
+ }
+ else
+ {
+ rightPositions.push_back(aPoint.x);
+ }
+ }
+ if(leftPositions.size() > 100)
+ {
+ isLeftNucleateExited = true;
+ }
+ if(rightPositions.size() > 100)
+ {
+ isRightNucleateExited = true;
+ }
+ if(isLeftNucleateExited && leftPositions.size() <= 80 &&
+ leftPositions.size() >= 10 && rightPositions.size() > 400 &&
+ leftPositions.size() > prevLeftSize)
+ {
+ int moleculeCnt(0);
+ for(unsigned int i(0); i != leftPositions.size(); ++i)
+ {
+ if(leftPositions[i] <= 0.5*theCenterPoint.x)
+ {
+ ++moleculeCnt;
+ }
+ }
+ cout << "moleculeCnt:" << moleculeCnt << endl;
+ int currStatus(4);
+ if(moleculeCnt < 0.5*leftPositions.size())
+ {
+ currStatus = 5; //Aberrant polar nucleation
+ }
+ if(isLeftNucleated)
+ {
+ Period = theTime - theLeftBeginTime;
+ ++theCycleCount;
+ theTotalPeriod += Period;
+ if(prevLeftStatus == 4 && currStatus == 4 && prevRightStatus == 4)
+ {
+ Status = 4; //Normal nucleations after the last complete cycle
+ }
+ else
+ {
+ Status = 5; //Aberrant polar nucleation
+ }
+ }
+ theLeftBeginTime = theTime;
+ isLeftNucleated = true;
+ isLeftNucleateExited = false;
+ prevLeftStatus = currStatus;
+ }
+ else if(isRightNucleateExited && rightPositions.size() <= 80 &&
+ rightPositions.size() >= 10 && leftPositions.size() > 400 &&
+ rightPositions.size() > prevRightSize)
+ {
+ int moleculeCnt(0);
+ for(unsigned int i(0); i != rightPositions.size(); ++i)
+ {
+ if(rightPositions[i] >= (theCenterPoint.x+0.5*theCenterPoint.x))
+ {
+ ++moleculeCnt;
+ }
+ }
+ cout << "moleculeCnt:" << moleculeCnt << endl;
+ int currStatus(4);
+ if(moleculeCnt < 0.5*rightPositions.size())
+ {
+ currStatus = 5; //Aberrant polar nucleation
+ }
+ if(isRightNucleated)
+ {
+ if(prevRightStatus == 4 && currStatus == 4 && prevLeftStatus == 4)
+ {
+ Status = 4; //Normal nucleations after the last complete cycle
+ }
+ else
+ {
+ Status = 5; //Aberrant polar nucleation
+ }
+ }
+ theRightBeginTime = theTime;
+ isRightNucleated = true;
+ isRightNucleateExited = false;
+ prevRightStatus = currStatus;
+ }
+ prevLeftSize = leftPositions.size();
+ prevRightSize = rightPositions.size();
+ cout << "theTime:" << theTime << " left:" << leftPositions.size() << " right:" << rightPositions.size() << " current period:" << Period << " avg period:" << theTotalPeriod/theCycleCount << endl;
+}
168 OscillationAnalysisProcess.hpp
@@ -0,0 +1,168 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __OscillationAnalysisProcess_hpp
+#define __OscillationAnalysisProcess_hpp
+
+#include <fstream> //provides ofstream
+#include <MethodProxy.hpp>
+#include "SpatiocyteProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_CLASS(OscillationAnalysisProcess, SpatiocyteProcess)
+{
+public:
+ LIBECS_DM_OBJECT(OscillationAnalysisProcess, Process)
+ {
+ INHERIT_PROPERTIES(Process);
+ PROPERTYSLOT_SET_GET(Real, MembraneBindTime);
+ PROPERTYSLOT_SET_GET(Real, LocalizationStartTime);
+ PROPERTYSLOT_SET_GET(Real, LocalizationEndTime);
+ PROPERTYSLOT_SET_GET(Real, OscillationTime);
+ PROPERTYSLOT_SET_GET(Real, AnalysisInterval);
+ PROPERTYSLOT_SET_GET(Real, Period);
+ PROPERTYSLOT_SET_GET(Integer, Status);
+ }
+ OscillationAnalysisProcess():
+ isLeftNucleated(false),
+ isLeftNucleateExited(true),
+ isRightNucleated(false),
+ isRightNucleateExited(true),
+ prevLeftSize(0),
+ prevRightSize(0),
+ Status(0),
+ theCycleCount(0),
+ AnalysisInterval(1),
+ LocalizationEndTime(13),
+ LocalizationEndTime2(5),
+ LocalizationStartTime(8),
+ LocalizationStartTime2(30),
+ MembraneBindTime(6),
+ OscillationTime(40) {}
+ virtual ~OscillationAnalysisProcess() {}
+ SIMPLE_SET_GET_METHOD(Real, MembraneBindTime);
+ SIMPLE_SET_GET_METHOD(Real, LocalizationStartTime);
+ SIMPLE_SET_GET_METHOD(Real, LocalizationEndTime);
+ SIMPLE_SET_GET_METHOD(Real, OscillationTime);
+ SIMPLE_SET_GET_METHOD(Real, AnalysisInterval);
+ SIMPLE_SET_GET_METHOD(Real, Period);
+ SIMPLE_SET_GET_METHOD(Integer, Status);
+ virtual void initializeLastOnce()
+ {
+ for(unsigned int i(0); i != theProcessSpecies.size(); ++i)
+ {
+ if(theProcessSpecies[i]->getVariable()->getID() == "MinD")
+ {
+ minD_m = theProcessSpecies[i];
+ }
+ else if(theProcessSpecies[i]->getVariable()->getID() == "MinDatp")
+ {
+ minD = theProcessSpecies[i];
+ }
+ else if(theProcessSpecies[i]->getVariable()->getID() == "MinEE")
+ {
+ minE = theProcessSpecies[i];
+ }
+ }
+ theCenterPoint = theSpatiocyteStepper->getCenterPoint();
+ }
+ virtual GET_METHOD(Real, StepInterval)
+ {
+ return MembraneBindTime;
+ }
+ virtual void fire()
+ {
+ if(theTime <= MembraneBindTime)
+ {
+ testMembraneBinding();
+ theTime = LocalizationStartTime;
+ }
+ else if(theTime >= LocalizationStartTime &&
+ theTime <= LocalizationEndTime )
+ {
+ testLocalization(2);
+ theTime += AnalysisInterval;
+ if(theTime > LocalizationEndTime)
+ {
+ theTime = LocalizationStartTime2;
+ }
+ }
+ else if(theTime >= LocalizationStartTime2 &&
+ theTime <= LocalizationEndTime2 )
+ {
+ testLocalization(3);
+ theTime += AnalysisInterval;
+ if(theTime > LocalizationEndTime2)
+ {
+ theTime = OscillationTime;
+ }
+ }
+ else
+ {
+ testOscillation();
+ theTime += AnalysisInterval;
+ }
+ thePriorityQueue->moveTop();
+ }
+protected:
+ virtual void testMembraneBinding();
+ virtual void testLocalization(int);
+ virtual void testOscillation();
+protected:
+ bool isLeftNucleated;
+ bool isLeftNucleateExited;
+ bool isRightNucleated;
+ bool isRightNucleateExited;
+ unsigned int prevLeftSize;
+ unsigned int prevRightSize;
+ int currStatus;
+ int prevLeftStatus;
+ int prevRightStatus;
+ int Status;
+ int theCycleCount;
+ double AnalysisInterval;
+ double LocalizationEndTime;
+ double LocalizationEndTime2;
+ double LocalizationStartTime;
+ double LocalizationStartTime2;
+ double MembraneBindTime;
+ double OscillationTime;
+ double Period;
+ double theLeftBeginTime;
+ double theRightBeginTime;
+ double theTotalPeriod;
+ Species* minD;
+ Species* minD_m;
+ Species* minE;
+ Point theCenterPoint;
+};
+
+#endif /* __OscillationAnalysisProcess_hpp */
33 PeriodicBoundaryDiffusionProcess.cpp
@@ -0,0 +1,33 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "PeriodicBoundaryDiffusionProcess.hpp"
+
+LIBECS_DM_INIT(PeriodicBoundaryDiffusionProcess, Process);
71 PeriodicBoundaryDiffusionProcess.hpp
@@ -0,0 +1,71 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __PeriodicBoundaryDiffusionProcess_hpp
+#define __PeriodicBoundaryDiffusionProcess_hpp
+
+#include "DiffusionProcess.hpp"
+
+LIBECS_DM_CLASS(PeriodicBoundaryDiffusionProcess, DiffusionProcess)
+{
+public:
+ LIBECS_DM_OBJECT(PeriodicBoundaryDiffusionProcess, Process)
+ {
+ INHERIT_PROPERTIES(DiffusionProcess);
+ }
+ PeriodicBoundaryDiffusionProcess() {}
+ virtual ~PeriodicBoundaryDiffusionProcess() {}
+ virtual void initialize()
+ {
+ DiffusionProcess::initialize();
+ theSpatiocyteStepper->setPeriodicEdge();
+ }
+ void initializeThird()
+ {
+ DiffusionProcess::initializeThird();
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->initMoleculeOrigins();
+ }
+ }
+ virtual void fire()
+ {
+ DiffusionProcess::fire();
+ for(vector<Species*>::const_iterator i(theProcessSpecies.begin());
+ i != theProcessSpecies.end(); ++i)
+ {
+ (*i)->relocateBoundaryMolecules();
+ }
+ }
+};
+
+#endif /* __PeriodicBoundaryDiffusionProcess_hpp */
257 PolymerFragmentationProcess.cpp
@@ -0,0 +1,257 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "PolymerFragmentationProcess.hpp"
+#include "SpatiocyteSpecies.hpp"
+
+LIBECS_DM_INIT(PolymerFragmentationProcess, Process);
+
+void PolymerFragmentationProcess::fire()
+{
+ //polyA + polyB -> (poly)C + (poly)D
+ if(A && B && C && D)
+ {
+ //A must be the reference subunit:
+ int reactantIndex(gsl_rng_uniform_int(getStepper()->getRng(),
+ theReactantSize));
+ Voxel* moleculeA(theReactants[reactantIndex]);
+ Voxel* moleculeC;
+ if(A->getVacantID() != C->getVacantID())
+ {
+ moleculeC = C->getRandomAdjoiningVoxel(moleculeA);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of A which can be occupied by C:
+ if(moleculeC == NULL)
+ {
+ return;
+ }
+ }
+ else
+ {
+ moleculeC = moleculeA;
+ }
+ Voxel* moleculeB(moleculeA->subunit->targetVoxels[theBendIndexA]);
+ Voxel* moleculeD;
+ if(B->getVacantID() != D->getVacantID())
+ {
+ //Find an adjoining vacant voxel of B belonging to the
+ //compartment of D which is not molecule C:
+ moleculeD = D->getRandomAdjoiningVoxel(moleculeB, moleculeC);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of B which can be occupied by D and the vacant voxel is not
+ //used by moleculeC:
+ if(moleculeD == NULL)
+ {
+ return;
+ }
+ }
+ else
+ {
+ moleculeD = moleculeB;
+ }
+ //When we remove the molecule A, this process' queue in priority
+ //queue will be updated since the substrate number has changed:
+ A->removeMolecule(moleculeA);
+ C->addMolecule(moleculeC);
+ B->removeMolecule(moleculeB);
+ D->addMolecule(moleculeD);
+ //Change the connections only after updating the queues of
+ //other depolymerize processes since the connections are used by
+ //the processes when checking:
+ if(C->getIsPolymer())
+ {
+ thePolymerizeProcess->removeContPoint(moleculeB->subunit,
+ &moleculeA->subunit->targetPoints[theBendIndexA]);
+ moleculeA->subunit->targetVoxels[theBendIndexA] = NULL;
+ }
+ else
+ {
+ thePolymerizeProcess->resetSubunit(moleculeA->subunit);
+ }
+ if(D->getIsPolymer())
+ {
+ moleculeB->subunit->sourceVoxels[theBendIndexB] = NULL;
+ }
+ else
+ {
+ thePolymerizeProcess->resetSubunit(moleculeB->subunit);
+ }
+ }
+ //polyA -> C + D
+ else if(A && !B && C && D)
+ {
+ //A must be the reference subunit:
+ int reactantIndex(gsl_rng_uniform_int(getStepper()->getRng(),
+ theReactantSize));
+ Voxel* moleculeA(theReactants[reactantIndex]);
+ Voxel* moleculeC;
+ Voxel* moleculeD;
+ if(A->getVacantID() != C->getVacantID())
+ {
+ moleculeC = C->getRandomAdjoiningVoxel(moleculeA);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of A which can be occupied by C:
+ if(moleculeC == NULL)
+ {
+ return;
+ }
+ if(A->getVacantID() != D->getVacantID())
+ {
+ //Find an adjoining vacant voxel of A belonging to the
+ //compartment of D which is not molecule C:
+ moleculeD = D->getRandomAdjoiningVoxel(moleculeA, moleculeC);
+ //Only proceed if we can find an adjoining vacant voxel
+ //of B which can be occupied by D and the vacant voxel is not
+ //used by moleculeC:
+ if(moleculeD == NULL)
+ {
+ moleculeD = D->getRandomAdjoiningVoxel(moleculeC);
+ if(moleculeD == NULL)
+ {
+ return;
+ }
+ }
+ }
+ else
+ {
+ moleculeD = moleculeA;
+ }
+ }
+ else
+ {
+ moleculeC = moleculeA;
+ moleculeD = D->getRandomAdjoiningVoxel(moleculeA);
+ if(moleculeD == NULL)
+ {
+ return;
+ }
+ }
+ //When we remove the molecule A, this process' queue in priority
+ //queue will be updated since the substrate number has changed:
+ A->removeMolecule(moleculeA);
+ C->addMolecule(moleculeC);
+ D->addMolecule(moleculeD);
+ //Change the connections only after updating the queues of
+ //other depolymerize processes since the connections are used by
+ //the processes when checking:
+ thePolymerizeProcess->resetSubunit(moleculeA->subunit);
+ }
+ for(vector<ReactionProcess*>::const_iterator
+ i(theInterruptingProcesses.begin());
+ i!=theInterruptingProcesses.end(); ++i)
+ {
+ (*i)->substrateValueChanged(theSpatiocyteStepper->getCurrentTime());
+ }
+}
+
+void PolymerFragmentationProcess::addSubstrateInterrupt(Species* aSpecies,
+ Voxel* aMolecule)
+{
+ //If first order:
+ if(!B)
+ {
+ addMoleculeA(aMolecule);
+ return;
+ }
+ //If second order:
+ else if(aSpecies == A)
+ {
+ Subunit* subunitA(aMolecule->subunit);
+ if(subunitA->targetVoxels.size() &&
+ subunitA->targetVoxels[theBendIndexA] &&
+ subunitA->targetVoxels[theBendIndexA]->id == B->getID())
+ {
+ addMoleculeA(aMolecule);
+ return;
+ }
+ }
+ if(aSpecies == B)
+ {
+ Subunit* subunitB(aMolecule->subunit);
+ if(subunitB->sourceVoxels.size() &&
+ subunitB->sourceVoxels[theBendIndexB] &&
+ subunitB->sourceVoxels[theBendIndexB]->id == A->getID())
+ {
+ addMoleculeA(subunitB->sourceVoxels[theBendIndexB]);
+ return;
+ }
+ }
+}
+
+void PolymerFragmentationProcess::removeSubstrateInterrupt(Species* aSpecies,
+ Voxel* aMolecule)
+{
+ if(aSpecies == A)
+ {
+ for(unsigned int i(0); i < theReactantSize; ++i)
+ {
+ if(aMolecule == theReactants[i])
+ {
+ theReactants[i] = theReactants[--theReactantSize];
+ substrateValueChanged(theSpatiocyteStepper->getCurrentTime());
+ return;
+ }
+ }
+ }
+ if(aSpecies == B)
+ {
+ Subunit* subunitB(aMolecule->subunit);
+ Voxel* moleculeA(subunitB->sourceVoxels[theBendIndexB]);
+ if(subunitB->sourceVoxels.size() && moleculeA &&
+ moleculeA->id == A->getID())
+ {
+ for(unsigned int i(0); i < theReactantSize; ++i)
+ {
+ if(moleculeA == theReactants[i])
+ {
+ theReactants[i] = theReactants[--theReactantSize];
+ substrateValueChanged(theSpatiocyteStepper->getCurrentTime());
+ return;
+ }
+ }
+ }
+ }
+}
+
+void PolymerFragmentationProcess::initializeLastOnce()
+{
+ A->addInterruptedProcess(reinterpret_cast<SpatiocyteProcess*>(this));
+ theBendIndexA = A->getBendIndex(BendAngle);
+ if(B)
+ {
+ if(A != B)
+ {
+ B->addInterruptedProcess(reinterpret_cast<SpatiocyteProcess*>(this));
+ }
+ theBendIndexB = B->getBendIndex(BendAngle);
+ }
+}
+
+
129 PolymerFragmentationProcess.hpp
@@ -0,0 +1,129 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// based on GillespieProcess.hpp
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+
+#ifndef __PolymerFragmentationProcess_hpp
+#define __PolymerFragmentationProcess_hpp
+
+#include <MethodProxy.hpp>
+#include "ReactionProcess.hpp"
+#include "PolymerizationProcess.hpp"
+
+class PolymerizationProcess;
+
+LIBECS_DM_CLASS(PolymerFragmentationProcess, ReactionProcess)
+{
+public:
+ LIBECS_DM_OBJECT(PolymerFragmentationProcess, Process)
+ {
+ INHERIT_PROPERTIES(ReactionProcess);
+ PROPERTYSLOT_SET_GET(Real, BendAngle);
+ }
+ SIMPLE_SET_GET_METHOD(Real, BendAngle);
+ PolymerFragmentationProcess():
+ theReactantSize(0),
+ BendAngle(0) {}
+ virtual ~PolymerFragmentationProcess() {}
+ virtual void initialize()
+ {
+ ReactionProcess::initialize();
+ if(getOrder() != 1 && getOrder() != 2)
+ {
+ THROW_EXCEPTION(ValueError,
+ String(getPropertyInterface().getClassName()) +
+ "[" + getFullID().asString() +
+ "]: This PolymerFragmentationProcess requires " +
+ "two substrates.");
+ }
+ if(p == -1)
+ {
+ p = k;
+ }
+ else
+ {
+ k = p;
+ }
+ }
+ virtual void initializeLastOnce();
+ virtual bool isContinuous() const
+ {
+ return true;
+ }
+ virtual GET_METHOD(Real, StepInterval)
+ {
+ return getPropensity()*
+ (-log(gsl_rng_uniform_pos(getStepper()->getRng())));
+ }
+ virtual void fire();
+ virtual void addSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule);
+ virtual void removeSubstrateInterrupt(Species* aSpecies, Voxel* aMolecule);
+ void addMoleculeA(Voxel* aMolecule)
+ {
+ ++theReactantSize;
+ if(theReactantSize > theReactants.size())
+ {
+ theReactants.push_back(aMolecule);
+ }
+ else
+ {
+ theReactants[theReactantSize-1] = aMolecule;
+ }
+ substrateValueChanged(theSpatiocyteStepper->getCurrentTime());
+ }
+ void setPolymerizeProcess(PolymerizationProcess* aProcess)
+ {
+ thePolymerizeProcess = aProcess;
+ }
+protected:
+ const double getPropensity() const
+ {
+ if(theReactantSize > 0 && p > 0)
+ {
+ return 1/(p*theReactantSize);
+ }
+ else
+ {
+ return libecs::INF;
+ }
+ }
+protected:
+ unsigned int theReactantSize;
+ int theBendIndexA;
+ int theBendIndexB;
+ double BendAngle;
+ PolymerizationProcess* thePolymerizeProcess;
+ vector<Voxel*> theReactants;
+};
+
+
+#endif /* __PolymerFragmentationProcess_hpp */
+
+
33 PolymerizationParameterProcess.cpp
@@ -0,0 +1,33 @@
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+// This file is part of E-Cell Simulation Environment package
+//
+// Copyright (C) 2006-2009 Keio University
+//
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//
+//
+// E-Cell is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// E-Cell is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+// See the GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public
+// License along with E-Cell -- see the file COPYING.
+// If not, write to the Free Software Foundation, Inc.,
+// 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+//
+//END_HEADER
+//
+// written by Satya Arjunan <satya.arjunan@gmail.com>
+// E-Cell Project, Institute for Advanced Biosciences, Keio University.
+//
+
+#include "PolymerizationParameterProcess.hpp"
+
+LIBECS_DM_INIT(PolymerizationParameterProcess, Process);
103 PolymerizationParameterProcess.hpp
@@ -0,0 +1,103 @@