squirmchemistry.cpp

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

// SquirmChemistry.cpp

#include "stdafx.h"
#include "SquirmChemistry.h"
#include "math.h"

bool SquirmChemistry::tryProteinArithmeticReactions(
	SquirmCell *cell)
{
	// below we assume first_enzyme=7, may actually be 8 or 9 or ...

	// if atom x has state i>=7 and type a,b or c
	// and is bonded to an atom in state 1 then apply this reaction:
	// x(i)y1 -> x0 + y(j)
	// where j = 3(i-7) + value(x) + 7
	// and value(a)=0, value(b)=1, value(c)=2

	// at the same time, search for where to apply this reaction:
	// d(i)y1 -> d(i)+y7

	if(cell->GetState()>=first_enzyme && 
		(cell->GetType()=='a' ||
		cell->GetType()=='b' || cell->GetType()=='c' ||
		cell->GetType()=='d'))
	{
		// search for a bonded atom in state 1
		POSITION pos = cell->GetBondedCells().GetHeadPosition();
		while(pos)
		{
			SquirmCell *y = cell->GetBondedCells().GetNext(pos);
			if(y->GetState()==1)
			{
				if(cell->GetType()!='d')
				{
					// apply the protein arithmetic reaction
					int i,j,value_of_x;
					i = cell->GetState();
					switch(cell->GetType())
					{
					case 'a': value_of_x=0; break;
					case 'b': value_of_x=1; break;
					case 'c': value_of_x=2; break;
					}
					j = 3*(i-first_enzyme)+value_of_x+first_enzyme;
					cell->SetState(0);
					y->SetState(j);
					cell->Debond(y);
					return true;
				}
				else { // type=='d'
					// enzyme produced, break off and continue priming (or stop if 'f')
                    if(y->GetType()!='f')
					    y->SetState(first_enzyme);
					cell->Debond(y);
					return true;
				}
			}
		}
	}

	return false;
}


bool SquirmChemistry::tryEnzymeReactions(SquirmCell *cell,
	const SquirmCellList& unbonded_neighbours)
{
	// for cells of type 'd' in state i>=first_enzyme, apply the enzyme reaction:
	
	//    d(i)   x(g)          d(i)   x(j)
	//          /b1/     ->          /b2/       
	//         y(h)                 y(k)

	// where:
	// N is the number of states (==first_enzyme)
	const int N = this->first_enzyme;
	// M is the number of types (6 for now)
	const int M = 6;

	// i = 2(2(M(M(N(N(Ng+h)+j)+k)+x)+y)+b1)+b2+first_enzyme
	//   = 4MMNNNg + 4MMNNh + 4MMNj + 4MMk + 4Mx + 4y + 2b1 + b2 + f_e
	// for x and y: 'a'=0, 'b'=1, 'c'=2, 'd'=3, 'e'=4, 'f'=5
	// for b1 and b2: unbonded=0 bonded=1

    // eg. a8 + a0 -> a8-a7 can be represented by: (N=18,M=6)
    // i = 4*6*6*18*18*18*8 + 4*6*6*18*18*0 + 4*6*6*18*8 + 4*6*6*7 + 4*6*0 + 4*0 + 2*0 + 1 + N
    //  = 6718464 + 0 + 20736 + 1008 + 0 + 0 + 0 + 1 + 18
    //  = 6740209 + 18
    // 6740209 is 110200102211101 in base 3 and so can be encoded as: bbacaabaccbbbab

    // eg. e9 + a8 -> e10-a8 can be represented by:
    // i = 4*6*6*18*18*18*9 + 4*6*6*18*18*8 + 4*6*6*18*10 + 4*6*6*8 + 4*6*4 + 4*0 + 2*0 + 1 + N
    //   = 7958689 + N
    // 7958689 is 112222100021021 in base 3, hence bbccccbaaacbacb

    if(cell->GetType()!='d' || cell->GetState()<first_enzyme)
		return false;

	unsigned long i=cell->GetState();
	int g,h,j,k,xval,yval,b1,b2;
	unsigned long divisor;
	
	// decode the reaction from the value of i

	i-=first_enzyme;

	divisor = 4*M*M*N*N*N;
	g = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 4*M*M*N*N;
	h = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 4*M*M*N;
	j = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 4*M*M;
	k = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 4*M;
	xval = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 4;
	yval = (int)floor(i/(float)divisor);
	i = i % divisor;

	divisor = 2;
	b1 = (int)floor(i/(float)divisor);
	i = i % divisor;

	b2 = i;


	// convert values xval and yval to types x and y
	char x,y;
	{ const char conv[6]={'a','b','c','d','e','f'}; 
	x=conv[xval]; y=conv[yval]; }

	// now search for the matching atoms nearby
	SquirmCellList x_candidates;
	getThoseOfTypeAndState(unbonded_neighbours,x,g,x_candidates);
	if(x_candidates.GetCount()<=0) return false;
	SquirmCellList y_candidates;
	getThoseOfTypeAndState(unbonded_neighbours,y,h,y_candidates);
	if(y_candidates.GetCount()<=0) return false;

	// search for a suitable pair
	POSITION pos = x_candidates.GetHeadPosition(),pos2;
	SquirmCell *x_candidate,*y_candidate;
	while(pos)
	{
		x_candidate = x_candidates.GetNext(pos);
		pos2 = y_candidates.GetHeadPosition();
		while(pos2)
		{
			y_candidate = y_candidates.GetNext(pos2);
			// is this a suitable pair?
			if(x_candidate != y_candidate && (
				(b1==0 && !x_candidate->HasBondWith(y_candidate)) ||
				(b1==1 && x_candidate->HasBondWith(y_candidate))))
			{
				// yes, apply the reaction
				x_candidate->SetState(j);
				y_candidate->SetState(k);
				if(b1==0 && b2==1)
					x_candidate->BondTo(y_candidate);
				else if(b1==1 && b2==0)
					x_candidate->Debond(y_candidate);
				return true;
			}
		}
	}
	
	return false;
}