squirmgrid.cpp

本程序模拟细胞的自我繁殖

// SquirmGrid.cpp

#include "stdafx.h"
#include "SquirmGrid.h"
#include "SquirmError.h"
#include "math.h"

#include "MainFrm.h"

//------------------------------------------------------------------------
SquirmGrid::SquirmGrid()
{
    this->created=false;

    // -- create the reactions --

	if(false) {
		// lipid chemistry
		//chemistry.addReaction(SquirmReaction('a',0,false,'b',0,1,true,1));

        chemistry.addReaction(SquirmReaction('a',0,false,'b',0,false,'d',1,false,1,true,1,false,1,false));
       
        //chemistry.addReaction(SquirmReaction('a',0,false,'b',0,false,'c',1,false,1,true,1,false,1,false));
	}

	if(true)
	{
		chemistry.SetFirstEnzymeState(18);
		
		// -- the reactions (2/3 nhood) --

        // attach the top to the membrane and trigger duplication
		chemistry.addReaction(SquirmReaction('e',9,false,'a',8,10,true,8)); 
		chemistry.addReaction(SquirmReaction('e',10,false,'e',0,false,'a',8,true,4,true,3,true,8,true));

		// duplication
		chemistry.addReaction(SquirmReaction('x',3,true,'x',4,true,'y',1,false,4,true,4,true,2,true));
        chemistry.addReaction(SquirmReaction('x',2,false,'x',0,3,true,5));
        chemistry.addReaction(SquirmReaction('x',3,true,'x',5,false,'y',4,true,4,true,3,true,4,false));

        // attach the bottom to the membrane and trigger splitting
		chemistry.addReaction(SquirmReaction('f',4,true,'f',3,false,'a',8,false,7,true,8,true,10,true));
		chemistry.addReaction(SquirmReaction('a',10,true,'a',8,false,'f',8,true,10,true,11,false,7,false));

		// splitting
        chemistry.addReaction(SquirmReaction('x',6,true,'y',4,1,true,7));
		chemistry.addReaction(SquirmReaction('x',7,true,'y',7,6,false,6));

        // join the membranes and release the information-carrying molecules
        chemistry.addReaction(SquirmReaction('e',6,true,'a',8,1,true,8));
		chemistry.addReaction(SquirmReaction('a',8,true,'e',1,true,'a',12,false,9,false,1,false,9,true));
		chemistry.addReaction(SquirmReaction('a',9,true,'a',9,true,'a',10,false,10,true,9,false,14,true));
		chemistry.addReaction(SquirmReaction('a',14,true,'a',10,true,'a',8,false,15,false,14,false,8,true));
		chemistry.addReaction(SquirmReaction('a',14,true,'e',1,false,'a',9,true,8,false,16,false,8,true));
        chemistry.addReaction(SquirmReaction('a',15,true,'e',1,8,false,16));

        // touch top to bottom to trigger protein-production (this is still necessary to change the f end)
        chemistry.addReaction(SquirmReaction('e',16,false,'f',1,17,false,17));

		// protein read-off (uses same 'zippering' mechanism as replication)
        chemistry.addReaction(SquirmReaction('e',17,true,'x',1,17,true,2));
        chemistry.addReaction(SquirmReaction('e',17,true,'x',3,true,'y',5,false,17,true,4,true,3,true));
        chemistry.addReaction(SquirmReaction('f',17,true,'x',4,true,'y',3,false,1,true,6,false,6,true));
		chemistry.addReaction(SquirmReaction('x',6,true,'e',17,true,'x',6,false,chemistry.GetFirstEnzymeState(),false,9,true,1,false));

		// -- very occasional mutations --

		if(false) {
			// occasional spontaneous formation of e1f1 pairs!
			chemistry.addReaction(SquirmReaction('e',0,false,'f',0,
				1,true,1,50000));
		}

		if(true) {
			// gene insertion happens sometimes (with the free atom states: 0)
			chemistry.addReaction(SquirmReaction('x',1,true,'y',1,false,'z',0,false,
				1,false,1,true,1,true,10000000));

			// gene deletion
			chemistry.addReaction(SquirmReaction('x',1,true,'y',1,true,'z',1,false,
				1,false,0,false,1,true,10000000));
		}
	}

	if(true)
	{
		// primitive motility reactions (just for state 1 right now)
		const int S=1;
		chemistry.addReaction(SquirmReaction('a',10,true,'a',11,false,'x',S,true,10,true,12,true,S,true));
		chemistry.addReaction(SquirmReaction('x',S,true,'y',S,false,'a',12,true,S,true,S,true,13,false));
		chemistry.addReaction(SquirmReaction('x',S,true,'a',13,true,'a',10,false,S,false,11,true,13,true));
		chemistry.addReaction(SquirmReaction('x',S,true,'y',S,true,'a',13,true,S,true,S,true,10,false));
	}

    if(true)
    {
        // a7 spikes get incorporated into an a8 membrane
        chemistry.addReaction(SquirmReaction('a',8,true,'a',7,false,'a',8,true,8,true,8,true,8,false));
    }

    if(false)
    {
        // a little experiment... not good with self-rep cell stuff but really nice on its own
        chemistry.addReaction(SquirmReaction('a',8,true,'a',8,true,'a',8,false,8,false,8,true,8,true));
    }

	if(false)
	{
		// membrane experiments 

		// catalyse a new membrane atom
		chemistry.addReaction(SquirmReaction('d',1,false,'a',0,1,false,1));
		// membrane growth
		chemistry.addReaction(SquirmReaction('a',2,true,'a',2,false,'a',1,false,2,false,2,true,2,true));
		// membrane merge
		chemistry.addReaction(SquirmReaction('d',2,false,'b',0,2,false,1));
		chemistry.addReaction(SquirmReaction('b',1,false,'a',2,true,'a',2,false,2,true,2,false,2,true));
		chemistry.addReaction(SquirmReaction('b',2,true,'a',2,false,'b',2,false,4,false,3,true,3,true));
		chemistry.addReaction(SquirmReaction('b',4,true,'b',3,true,'a',2,false,5,false,5,false,3,true));
		//chemistry.addReaction(SquirmReaction('b',5,true,'a',3,2,true,2));
		chemistry.addReaction(SquirmReaction('a',2,true,'b',5,true,'a',3,false,2,false,0,false,2,true));
	}

    srand(time(NULL));

    this->do_flood = true;
    this->FLOOD_PERIOD = 50000;
    this->flood_sector = 2; // which quarter to flood next

	this->do_swap=false; // legacy code from ECAL03 experiments
	this->SWAP_PERIOD=this->FLOOD_PERIOD*5;

	this->do_cosmic=false; // cosmic-ray events flip the state of an atom at random
	this->cosmic_cases = 100000;

	this->do_output=false; // unless can open outputs
    if(false) // set to true and enter a filename below to get a log of reactions
              // (see RecomputeVelocitiesAndReact() )
    {
        this->output = fopen("C:\\Tim\\Papers\\Cells2004\\output\\log.csv","wt");
        if(this->output)
        {
            this->do_output = true;
            fprintf(this->output,"iterations,reaction,molecule\n");
        }
    }

    this->probe = NULL;

    this->OFFSET = 3;
}

/* old reaction scheme (neat, works well in continuous space)

  		// duplication
		chemistry.addReaction(SquirmReaction('e',1,false,'e',0,4,true,3));
		chemistry.addReaction(SquirmReaction('x',3,true,'x',4,true,'y',1,false,4,true,4,true,2,true));
        chemistry.addReaction(SquirmReaction('x',2,false,'x',0,3,true,5));
        chemistry.addReaction(SquirmReaction('x',3,true,'x',5,false,'y',4,true,4,true,3,true,4,false));

		// splitting
		chemistry.addReaction(SquirmReaction('f',3,true,'f',4,6,false,6));
		chemistry.addReaction(SquirmReaction('x',6,true,'y',4,1,true,7));
		chemistry.addReaction(SquirmReaction('x',7,true,'y',7,6,false,6));

		// protein read-off (uses same 'zippering' mechanism as replication)
		chemistry.addReaction(SquirmReaction('e',6,true,'x',1,false,'x',0,false,5,true,4,true,3,true));
		chemistry.addReaction(SquirmReaction('x',6,true,'e',5,true,'x',6,false,chemistry.GetFirstEnzymeState(),false,1,true,1,false));

*/

void SquirmGrid::Create(float x,float y)
{
    if(this->created)
    {
        this->~SquirmGrid();
    }

	this->size.x = x;
	this->size.y = y;

	// how many slots do we need?
	this->slot_size = RADIUS*5.0F;
	this->slots_X = (int)ceil(this->size.x / this->slot_size);
	this->slots_Y = (int)ceil(this->size.y / this->slot_size);

    // create a 2D array of SquirmCellSlot's NR-stylee
    SquirmCellSlot *matrix = new SquirmCellSlot[this->slots_X*this->slots_Y];
    this->cell_grid = new SquirmCellSlot*[this->slots_X];
    for(int i = 0; i < this->slots_X; i++)
        this->cell_grid[i] = &matrix[i*this->slots_Y];

    this->created=true;

    this->iterations=0;

}

SquirmCellSlot* SquirmGrid::GetSlot(C2dVector p)
{
	int x = (int)floor(p.x/slot_size);
	int y = (int)floor(p.y/slot_size);
	if(x>=0 && y>=0 && x<this->slots_X && y<this->slots_Y)
		return &this->cell_grid[x][y];
	else
	{
		AfxMessageBox("Error!");
		ASSERT(false);
		return NULL;
	}
}

//------------------------------------------------------------------------
SquirmGrid::~SquirmGrid()
{
  delete [] *this->cell_grid;
  delete [] this->cell_grid;

  delete this->probe;

  // delete all the cells
  POSITION pos = this->cell_list.GetHeadPosition();
  while(pos)
  {
      delete this->cell_list.GetNext(pos);
  }

  this->created=false;

  if(this->do_output)
    fclose(this->output);
}


void SquirmGrid::InitSimple(int n)
{
    // initialise a replicator
    if(true)
    {
        //bbacaabaccbbbab, bbccccbaaacbacb
        CreateNewCell(20.0F+1*RADIUS*2.0F,20.0F,'e',9)->BondTo(
        CreateNewCell(20.0F+2*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+3*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+4*RADIUS*2.0F,20.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+5*RADIUS*2.0F,20.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+6*RADIUS*2.0F,20.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+7*RADIUS*2.0F,20.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+8*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+9*RADIUS*2.0F,20.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+10*RADIUS*2.0F,20.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+11*RADIUS*2.0F,20.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+12*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+13*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+14*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+15*RADIUS*2.0F,20.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+16*RADIUS*2.0F,20.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+17*RADIUS*2.0F,20.0F,'d',1)->BondTo(
        CreateNewCell(20.0F+18*RADIUS*2.0F,20.0F,'f',1))))))))))))))))));
    }

    if(false)
    {
        CreateNewCell(20.0F+1*RADIUS*2.0F,160.0F,'c',9)->BondTo(
        CreateNewCell(20.0F+2*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+3*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+4*RADIUS*2.0F,160.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+5*RADIUS*2.0F,160.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+6*RADIUS*2.0F,160.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+7*RADIUS*2.0F,160.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+8*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+9*RADIUS*2.0F,160.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+10*RADIUS*2.0F,160.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+11*RADIUS*2.0F,160.0F,'c',1)->BondTo(
        CreateNewCell(20.0F+12*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+13*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+14*RADIUS*2.0F,160.0F,'b',1)->BondTo(
        CreateNewCell(20.0F+15*RADIUS*2.0F,160.0F,'a',1)->BondTo(
        CreateNewCell(20.0F+16*RADIUS*2.0F,160.0F,'d',1)->BondTo(
        CreateNewCell(20.0F+18*RADIUS*2.0F,160.0F,'f',1)))))))))))))))));
    }

    // create a catalyst
    if(false)
    {
        CreateNewCell(10,10,'d',1);
    }

    	// initialise a loop of a8's
    if(true)
	{
		int length=16;
		int width=25;
		float start_x=RADIUS,start_y=RADIUS;
		const int S=8;
		SquirmCell *start =  CreateNewCell(start_x,start_y,'a',S);
		SquirmCell *a=start,*b;
        int i;
		// a string down the left
		for(i=0;i<length;i++)
		{
			b = CreateNewCell(start_x,start_y+(i+1)*RADIUS*2.0F,'a',S);
			b->BondTo(a);
			a=b;
		}
		// along the bottom
		for(i=0;i<width;i++)
		{
			b = CreateNewCell(start_x+(i+1)*RADIUS*2.0F,start_y+length*RADIUS*2.0F,'a',S);
			b->BondTo(a);