Team:Osaka/Modeling

From 2013.igem.org

Revision as of 03:56, 28 September 2013 by Kakimoto (Talk | contribs)

team team team team team team team team
team
iGEM Osaka Theory Group have made mathematical models to understand how colony patterns look like.

Syuichi Iwakiri built differential equations to describe the movement of E.coli and nutrition.

And Hiroki Nishiyama did simulations of the equations.

In the equation on the picture , A means the number of type A-Ecoli , B means the number of type B-Ecoli ,

and N means the concentration of nutrition.






The first partial differential equation of A has 3 factors.

The first one is the diffusion of A , the second one is the growth of A ,

and the third one is the transition from A to B. The transition happens only when N is smaller than No.

The second partial differential equation of B has 2 factors.

The first one is the diffusion of B , and the second one is the transition from A to B.

The third partial differential equation of N has 3 factors.

The first one is the diffusion of N , and the second one is the consume by A ,

the third one is the virtual consume by B (b is the consume rate and c is the produce rate).

The Theory Group have discovered a unique patterns ,

for example , a fractal pattern and cross-like pattern.



We show several patterns that we got from the simulations below.

The gif files show how the number of A type bacterias(Z-axis) develop in 2D area(X-Y plate).


module my_precision

 implicit none
   integer,parameter :: wp = selected_real_kind(p=20)  !Change P to change max digit
   

end module my_precision

program abc2

use my_precision
implicit none

 real(wp)::A(200,200,2000),B(200,200,2000),N(200,200,2000)

!$$$$$$  !RESULT(2000,2000,2000)

REAL::Da,Db,Dn,N1,N2,N3,Ma,Mb,Ea,Eb,F,Ra,Rb,K,C
INTEGER::T,X,Y,Q,R,S,dT,dX,dY,TEND,XEND,YEND,PT

print*, 'start'

!............INPUT_ZERO.............. T=0 X=0 Y=0 dT=0 dX=0 dY=0 Q=0 R=0 S=0 TEND=0 XEND=0 YEND=0 Da=0 Db=0 Dn=0 PT=0.d0 N1=0.d0 N2=0.d0 N3=0.d0 Ma=0.d0 Mb=0.d0 Ea=0.d0 Eb=0.d0 F=0.d0 Ra=0.d0 Rb=0.d0 C=0.d0 K=0.d0


DO S=1,2000

 PRINT*,'S',S
 DO R=1,200
   print*,'R',R
   DO Q=1,200
       A(Q,R,S)=0 !A(X,Y,T)
       B(Q,R,S)=0
       N(Q,R,S)=0
       !RESULT(2000,2000,2000)=0
   ENDDO
 ENDDO

ENDDO

print*,'finish input zero'

!............INITIALIZATION..........

!..Mortaility rate.. Ma=0.05 !0.05 Mb=0.05 !0.05

!..Energy Usage Rate.. Ea=0.145 !0.145 Eb=0.145

!..Energy Produce Rate.. F=0.25 !0.25 0.5?

!..Growth Rate.. Ra=0.2 !0.2 Rb=0.2  !0.2

!..A to B Transformation Rate.. C=0.05 !0.05

!..Kankyou Shuyousu.. K=1000

!..AREA.. TEND=495 XEND=100 YEND=100

!..delta..do not make them 0.1 dX=1 dY=1 dT=1

!..Kakusan Keisuu.. Da=0.1 Db=0.1 Dn=0.1

A(50,50,3)=100 B(50,50,3)=100

 DO Q=3,XEND+1
   DO R=3,YEND+1
       N(Q,R,3)=100
   ENDDO
 ENDDO

N1=20 N2=60



!............MAIN CALCULATION........


DO T=3,TEND

   print*,'T',T,DX,DX    
 DO X=3,XEND
   DO Y=3,YEND
     
     PT=T+1
     
     IF (N(X,Y,T) < N1) THEN
        A(X,Y,PT)=Da*dT*((A(X+dX,Y,T)+A(X-dX,Y,T)-2*A(X,Y,T))/dX**2+(A(X,Y+dY,T)+A(X,Y-dY,T)-2*A(X,Y,T))/dY**2)+A(X,Y,T)
        A(X,Y,PT)=A(X,Y,PT)-Ma*A(x,y,t)-C*A(x,y,t)
        
        B(X,Y,PT)=Db*dT*((B(X+dX,Y,T)+B(X-dX,Y,T)-2*B(X,Y,T))/dX**2+(B(X,Y+dY,T)+B(X,Y-dY,T)-2*B(X,Y,T))/dY**2)+B(X,Y,T)
        B(X,Y,PT)=B(X,Y,PT)-Mb*B(x,y,t)+C*A(x,y,t)
        
        N(X,Y,PT)=Dn*dT*((N(X+dX,Y,T)+N(X-dX,Y,T)-2*N(X,Y,T))/dX**2+(N(X,Y+dY,T)+N(X,Y-dY,T)-2*N(X,Y,T))/dY**2)+N(X,Y,T)
        N(X,Y,PT)=N(X,Y,PT)-Ea*A(x,y,t)+(F-Eb)*B(x,y,t)
        
     ENDIF
     IF (N(X,Y,T) >= N1 .and. N(X,Y,T) <N2) THEN
        A(X,Y,PT)=Da*dT*((A(X+dX,Y,T)+A(X-dX,Y,T)-2*A(X,Y,T))/dX**2+(A(X,Y+dY,T)+A(X,Y-dY,T)-2*A(X,Y,T))/dY**2)+A(X,Y,T)
        A(X,Y,PT)=A(X,Y,PT)+Ra*A(x,y,t)*(1-A(x,y,t)/K)-C*A(x,y,t)
        
        B(X,Y,PT)=Db*dT*((B(X+dX,Y,T)+B(X-dX,Y,T)-2*B(X,Y,T))/dX**2+(B(X,Y+dY,T)+B(X,Y-dY,T)-2*B(X,Y,T))/dY**2)+B(X,Y,T)
        B(X,Y,PT)=B(X,Y,PT)+C*A(x,y,t)!+Rb*B(x,y,t)*(1-B(x,y,t)/K)
        
        N(X,Y,PT)=Dn*dT*((N(X+dX,Y,T)+N(X-dX,Y,T)-2*N(X,Y,T))/dX**2+(N(X,Y+dY,T)+N(X,Y-dY,T)-2*N(X,Y,T))/dY**2)+N(X,Y,T)
        N(X,Y,PT)=N(X,Y,PT)-Ea*A(x,y,t)+(F-Eb)*B(x,y,t)
        
     ENDIF
     IF (N(X,Y,T) >= N2) THEN
        A(X,Y,PT)=Da*dT*((A(X+dX,Y,T)+A(X-dX,Y,T)-2*A(X,Y,T))/dX**2+(A(X,Y+dY,T)+A(X,Y-dY,T)-2*A(X,Y,T))/dY**2)+A(X,Y,T)
        A(X,Y,PT)=A(X,Y,PT)+Ra*A(x,y,t)*(1-A(x,y,t)/K)
        
        B(X,Y,PT)=Db*dT*((B(X+dX,Y,T)+B(X-dX,Y,T)-2*B(X,Y,T))/dX**2+(B(X,Y+dY,T)+B(X,Y-dY,T)-2*B(X,Y,T))/dY**2)+B(X,Y,T)
        !B(X,Y,PT)=B(X,Y,PT)+Rb*B(x,y,t)*(1-B(x,y,t)/K)
        
        N(X,Y,PT)=Dn*dT*((N(X+dX,Y,T)+N(X-dX,Y,T)-2*N(X,Y,T))/dX**2+(N(X,Y+dY,T)+N(X,Y-dY,T)-2*N(X,Y,T))/dY**2)+N(X,Y,T)
        N(X,Y,PT)=N(X,Y,PT)-Ea*A(x,y,t)+(F-Eb)*B(x,y,t)
     ENDIF
     
   ENDDO
 ENDDO

ENDDO


OPEN(22,FILE='Blast21c.txt',POSITION='APPEND')

 DO X=3,XEND
   DO Y=3,YEND
     WRITE(22,*)X,Y,B(X,Y,905),B(X,Y,805),B(X,Y,705),B(X,Y,605),B(X,Y,505),B(X,Y,405),B(X,Y,305),B(X,Y,205),B(X,Y,105),B(X,Y,5)
     !B(X,Y,905),B(X,Y,895),B(X,Y,815),B(X,Y,785),B(X,Y,755),B(X,Y,725),B(X,Y,695),B(X,Y,665),B(X,Y,635),B(X,Y,605)
   ENDDO
 ENDDO

CLOSE(22)

OPEN(24,FILE='Alast21c.txt',POSITION='APPEND')

 DO X=3,XEND
   DO Y=3,YEND
     WRITE(24,*)X,Y,A(X,Y,905),A(X,Y,805),A(X,Y,705),A(X,Y,605),A(X,Y,505),A(X,Y,405),A(X,Y,305),A(X,Y,205),A(X,Y,105),A(X,Y,5)
     
     !X,Y,A(X,Y,905),A(X,Y,895),A(X,Y,815),A(X,Y,785),A(X,Y,755),A(X,Y,725),A(X,Y,695),A(X,Y,665),A(X,Y,635),A(X,Y,605)
   ENDDO
 ENDDO

CLOSE(24)


OPEN(26,FILE='Alast21d.txt',POSITION='APPEND')

 DO X=3,XEND
   DO Y=3,YEND

WRITE(26,*)X,Y,A(X,Y,1785),A(X,Y,1705),A(X,Y,1605),A(X,Y,1505),A(X,Y,1405),A(X,Y,1305),A(X,Y,1205),A(X,Y,1055),A(X,Y,1005)

   ENDDO
 ENDDO

CLOSE(26)

OPEN(28,FILE='Blast21d.txt',POSITION='APPEND')

 DO X=3,XEND
   DO Y=3,YEND

WRITE(28,*)X,Y,B(X,Y,1785),B(X,Y,1705),B(X,Y,1605),B(X,Y,1505),B(X,Y,1405),B(X,Y,1305),B(X,Y,1205),B(X,Y,1055),B(X,Y,1005)

   ENDDO
 ENDDO

CLOSE(28)

OPEN(32,FILE='Nlast21.txt',POSITION='APPEND')


 DO X=3,XEND
   DO Y=3,YEND
     WRITE(32,*)X,Y,N(X,Y,1995),N(X,Y,275),N(X,Y,255),N(X,Y,235),N(X,Y,215),N(X,Y,195),N(X,Y,175),N(X,Y,155),N(X,Y,135),N(X,Y,115)
     
   ENDDO
 ENDDO


CLOSE(32)


end program abc2