"
Team:Osaka/Modeling
From 2013.igem.org
Poribukuro (Talk | contribs) |
|||
| Line 87: | Line 87: | ||
</html> | </html> | ||
| - | + | BELOW WE SHOW ONE OF THE PROGRAMMING CODE WE HAVE MADE. <br><br> | |
module my_precision | module my_precision | ||
Revision as of 04:02, 28 September 2013
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).
BELOW WE SHOW ONE OF THE PROGRAMMING CODE WE HAVE MADE.
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
