Subroutine mc_k_ion > (XiK, YiK, VxK, VyK, XiNa, YiNa, VxNa, VyNa, > walltop, wallbottom, poreleft, poreright, nrIter, > K_flux, Na_flux, K_hits, K_bind, nr_ions, > P_open, ProbLeave) C C SUBOUTINE MC_K_ION C C (C) Antonius M.J. VanDongen, 2003 C C This subroutine implements the Monte Carlo simulation of the 1-site C model described in the accompanying paper. It is written in FORTRAN-77 C and was complied to produce a DLL that can be called from a Visual Basic C program that implemented the graphical user interface. C Before calling the routine, the X,Y coordinates and velocities must be C initialized to random values. C C ------------------------------------------------------------------------------ C INPUT PARAMETERS C ------------------------------------------------------------------------------ C XiK, YiK Arrays with X,Y Coordinates of K ions C VxK, VyK Arrays with velocites of K ions C XiNa, YiNa Arrays with X,Y Coordinates of Na ions C VxNa, VyNa Arrays with velocites of Na ions C C walltop, wallbottom Y-ccordinates of the compartment top and bottom C The left/right border of the compartment are 10 & 9990, resp C C poreleft, poreright Y-coordinates of the filter entrance C nrIter Number of iterations C nr_ions (0:2) nr_ions(1) = number of K ions, C nr_ions(2) = number of Na ions C P_open Probability of being open (low affinity) = Plow C ProbLeave Probability of leaving the filter (in the low affinity state) C C ------------------------------------------------------------------------------ C OUTPUT PARAMETERS C ------------------------------------------------------------------------------ C K_flux, Total number of K ions that permeated C Na_flux Total number of Na ions that permeated C K_hits Total number of K ions that hit the filter entrance C K_bind Total number of K ions that was trapped in the filter C C ------------------------------------------------------------------------------ C INTERNAL PARAMETERS C ------------------------------------------------------------------------------ C occupied =0 when filter is unoccupied C >0 when occupied with K ion (value is K ion number) C <0 when occupied with Na ion (value is -1 * Na ion number) C integer*4 NrIter integer*2 xik(0:*),yik(0:*), vxk(0:*), vyk(0:*) integer*2 xina(0:*),yina(0:*), vxna(0:*), vyna(0:*) integer*2 walltop, wallbottom, poreleft, poreright integer*2 K_flux, Na_flux, K_hits, K_bind integer*2 occupied, gate, nr_ions(0:2) real*4 P_open, ProbLeave call one_site > (XiK, YiK, VxK, VyK, XiNa, YiNa, VxNa, VyNa, > walltop, wallbottom, poreleft, poreright, nrIter, > K_flux, Na_flux, K_hits, K_bind, gate, nr_ions, > P_open, occupied, ProbLeave) return end subroutine one_site > (XiK, YiK, VxK, VyK, XiNa, YiNa, VxNa, VyNa, > walltop, wallbottom, poreleft, poreright, nrIter, > K_flux, Na_flux, K_hits, K_bind, gate, nr_ions, > P_open, occupied, ProbLeave) integer*4 NrIter, t integer*2 xik(0:*),yik(0:*), vxk(0:*), vyk(0:*) integer*2 xina(0:*),yina(0:*), vxna(0:*), vyna(0:*) integer*2 walltop, wallbottom, poreleft, poreright integer*2 K_flux, Na_flux, K_hits, K_bind integer*2 occupied, gate, nr_ions(0:2) real*4 P_open, po1, ProbLeave real*4 middle integer*2 thisx, thisy, nextx, nexty integer*2 origin1, destination integer*2 j, top integer*4 tsum integer*4 ksum c calcualte position of the membrane middle = wallbottom - 0.5*(abs(wallbottom-walltop)) Po1 = P_open tsum = 0 ksum = 0 idum = -1000 * rand() c Outer loop: Do nrIter time steps for K and Na ions do t = 1, nrIter if (ran2(idum) .gt. Po1) then c if random number(0:1) is larger then Plow then move to high affinty if (occupied.lt.0)then c if filter is occupied by a Na ion, it can not move to high affinity gate = 1 else gate = 0 endif else gate = 1 endif c K ion loop: update all K ions do j=1, nr_ions(1) c j = current K ion thisx = XiK(j) thisy = YiK(j) c (thisx, thisy) = current X,Y coordinates of K ion j nexty = YiK(j) + VyK(j) nextx = XiK(j) + VxK(j) c (nextx, nexty) = new X,Y coordinates of K ion j c check for collisions with outside box If (nextx.le.10.or.nextx.ge.9990.or. > nexty.le.10.or.nexty.ge.9990) Then c K ion leaves box, replace with new K ion If (nexty .lt. middle) Then top = 1 else top = 0 endif call new_K_ion + (j, XiK, YiK, VxK, VyK, wallbottom, walltop, top) goto 10 endif c First check whether ion stays in solution, if not, replace If (thisy.lt.walltop .and. nexty .lt. walltop) goto 10 If (thisy.gt.wallbottom .and. nexty .gt. wallbottom) goto 10 c If (thisy .lt. walltop .and. nexty .ge. walltop) then c K ion is in external solution and enters membrane zone If (nextx .ge. PoreLeft .and. nextx .le. PoreRight) Then c K ion enters pore region K_hits = K_hits + 1 If (occupied .eq. 0) Then c site 1 is empty: trap K ion VxK(j) = 0 VyK(j) = 0 XiK(j) = poreleft + 0.44 * abs(poreright-poreleft) YiK(j) = walltop + 0.40 * abs(wallbottom-walltop) occupied = j K_bind = K_bind + 1 origin1 = +1 goto 10 End If Endif C Site is occupied or K ion missed the pore: bounce current K ion VyK(j) = -VyK(j) goto 10 Endif if (thisy .gt. wallbottom .and. nexty .le. wallbottom) then c K ion is in internal solution and enters membrane zone If (nextx .ge. PoreLeft .and. nextx .le. PoreRight) Then c K ion will enter pore region If (occupied .eq. 0) Then c site 1 is empty VxK(j) = 0 VyK(j) = 0 XiK(j) = poreleft + 0.44 * abs(poreright-poreleft) YiK(j) = walltop+0.40 * abs(wallbottom-walltop) occupied = j origin1 = -1 goto 10 endif endif c Site is occupied or K ion missed the pore: bounce current K ion VyK(j) = -VyK(j) goto 10 endif If (occupied .eq. j)then c Current K ion is in site 1, check whether it is moving if (gate .eq. 1) then prob_leave = ProbLeave else goto 10 endif If (ran2(idum).lt. prob_leave) then c K ion is leaving site 1 VxK(j) = 10+400 * (ran2(idum) - 0.5) if (ran2(idum).lt. 0.5) then VyK(j) = -10 - 200 * ran2(idum) destination = +1 else VyK(j) = 10 + 200 * ran2(idum) destination = -1 endif occupied = 0 c If K ion came from internal solution (origin=-1), then count as k flux if ((origin1 .eq. -1 .and. destination.eq. +1) .or. + (origin1 .eq. +1 .and. destination.eq. -1)) then K_flux=K_flux+1 ksum = ksum + 1 if (destination.eq.1) then top = 0 else top = 1 endif call new_K_ion + (j, XiK, YiK, VxK, VyK, wallbottom, walltop, top) goto 10 endif origin1 = 0 goto 10 Endif Endif 10 continue end do c Na ion loop: update all Na ions do j=1, nr_ions(2) thisx = XiNa(j) thisy = YiNa(j) nexty = YiNa(j) + VyNa(j) nextx = XiNa(j) + VxNa(j) c check for collisions with outside box if (nextx.le.10.or.nextx.ge.9990.or. > nexty.le.10.or.nexty.ge.9990) then c Na ion leaves box, replace with new Na ion if (nexty .lt. middle) then top = 1 else top = 0 endif call new_Na_ion + (j, XiNa, YiNa, VxNa, VyNa, wallbottom, walltop, top) goto 20 endif c First check whether ion stays in solution If (thisy.lt.walltop .and. nexty .lt. walltop) goto 20 If (thisy.gt.wallbottom .and. nexty .gt. wallbottom) goto 20 c if (thisy .lt. walltop .and. nexty .ge. walltop) then c Na ion is in external solution and enters membrane zone If (nextx .ge. PoreLeft .and. nextx .le. PoreRight) Then c Na ion enters pore region Na_hits = Na_hits + 1 If (occupied .eq. 0) Then c site 1 is empty If (gate.eq.0) then c Filter in high-affinity state: bounce Na ion VyNa(j) = -VyNa(j) Else c filter in low affinity state: Na ion enters filter VxNa(j) = 0 VyNa(j) = 0 XiNa(j) = poreleft + 0.44*abs(poreright-poreleft) YiNa(j) = walltop + 0.40*abs(wallbottom-walltop) occupied = -j Na_bind = Na_bind+1 origin1 = +1 Endif goto 20 End If Endif C Site is occupied or Na ion missed the pore: bounce current Na ion VyNa(j) = -VyNa(j) goto 20 endif if (thisy .gt. wallbottom .and. nexty .le. wallbottom) then c Na ion is in internal solution and enters membrane zone If (nextx .ge. PoreLeft .and. nextx .le. PoreRight) Then c Na ion will enter pore region If (occupied .eq. 0) Then c site 1 is empty If (gate.eq.0) then VyNa(j) = -VyNa(j) Else VxNa(j) = 0 VyNa(j) = 0 XiNa(j) = poreleft + 0.44*abs(poreright-poreleft) YiNa(j) = walltop+0.40*abs(wallbottom-walltop) occupied = -j origin1 = -1 endif goto 20 endif endif c Site is occupied or Na ion missed the pore: bounce current Na ion VyNa(j) = -VyNa(j) goto 20 endif If (occupied .eq. -j)then c Current Na ion is in site 1, check whether it is moving if (gate .eq. 1) then prob_leave = ProbLeave else goto 20 endif If (ran2(idum).lt. prob_leave) then c Na ion is leaving site 1 c Na ion moves into internal or external solution VxNa(j) = 10+400 * (ran2(idum) - 0.5) if (ran2(idum) .lt. 0.5) then VyNa(j) = -10 - 200 * ran2(idum) destination = +1 else VyNa(j) = 10 + 200 * ran2(idum) destination = -1 endif occupied = 0 c If Na ion came from internal solution (origin=-1), then count as Na flux If ((origin1 .eq. -1 .and. destination.eq. +1) .or. + (origin1 .eq. +1 .and. destination.eq. -1) ) then Na_flux=Na_flux+1 If (destination.eq.1) then top = 0 Else top = 1 Endif call new_Na_ion + (j, XiNa, YiNa, VxNa, VyNa, wallbottom, walltop, top) goto 20 endif origin1 = 0 goto 20 Endif Endif 20 continue end do 30 continue do j=1,nr_ions(1) XiK(j)=XiK(j)+VxK(j) YiK(j)=YiK(j)+VyK(j) end do do j=1,nr_ions(2) XiNa(j)=XiNa(j)+VxNa(j) YiNa(j)=YiNa(j)+VyNa(j) enddo end do return end FUNCTION ran2(idum) INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END subroutine new_K_ion + (j, XiK, YiK, VxK, VyK, wallbottom, walltop, top) integer*2 XiK(0:*), YiK(0:*), VxK(0:*), VyK(0:*) integer*2 walltop, wallbottom integer*2 j, top if (top .eq. 1) then XiK(j) = 10+ran2(idum)*9980 YiK(j) = ran2(idum)*0.90*walltop VxK(j) = 1+400 * (ran2(idum) - 0.5) VyK(j) = 1+400 * (ran2(idum) - 0.5) else XiK(j) = 10+ran2(idum)*9980 YiK(j) = wallbottom+500+ran2(idum)*(9500-wallbottom) VxK(j) = 1+400 * (ran2(idum) - 0.5) VyK(j) = 1+400 * (ran2(idum) - 0.5) endif c release "frozen" ions if (VxK(j) .eq. 0 .and. VyK(j) .eq. 0) then VxK(j) = 1 VyK(j) = 1 endif return end subroutine new_Na_ion + (j, XiNa, YiNa, VxNa, VyNa, wallbottom, walltop, top) integer*2 XiNa(0:*), YiNa(0:*), VxNa(0:*), VyNa(0:*) integer*2 walltop, wallbottom integer*2 j, top if (top .eq. 1) then XiNa(j) = 10+ran2(idum)*9980 YiNa(j) = ran2(idum)*0.90*walltop VxNa(j) = 1+400 * (ran2(idum) - 0.5) VyNa(j) = 1+400 * (ran2(idum) - 0.5) else XiNa(j) = 10+ran2(idum)*9980 YiNa(j) = wallbottom+500+ran2(idum)*(9500-wallbottom) VxNa(j) = 1+400 * (ran2(idum) - 0.5) VyNa(j) = 1+400 * (ran2(idum) - 0.5) endif c release "frozen" ions if (VxNa(j) .eq. 0 .and. VyNa(j) .eq. 0) then VxNa(j) = 1 VyNa(j) = 1 endif return end