Team:Osaka/Modeling

From 2013.igem.org

(Difference between revisions)
Line 94: Line 94:
</body>
</body>
</html>
</html>
 +
 +
 +
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

Revision as of 03:56, 28 September 2013

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

Retrieved from "http://2013.igem.org/Team:Osaka/Modeling"