% % File ~/CellExperiments/Dextran/Phalloidin/0422Movingf3.jpg % 1 pixel = 1/6 micrometer % L scale = 10 micrometer % t scale = 40 sec % V scale = 0.25 micrometer/sec % % Solution of the nonlinear problem % % \ksi0 F^2 (Vf - Vp) + \nabla p = 0, (1) % \nabla \cdot (h*[Vf*ph+(1-ph)*Vp]) = 0, (2) % Fs - \nabla \cdot (Vp*F) - \nabla \cdot (Vf*G) +D\nabla^2*G = 0 (3) % % % where \nabla={d/dx, d/dy}, % Vp = {Vpx, Vpy} is the F-actin velocity field (external data). % Vf={Vfx, Vfy} is 2D vector fluid velocity field, % F is filament(F)-actin density (external data), Fs - is its uniform source % p is fluid pressure, % G is density of G-actin % h(x,y) is the height of the cell, % ph=1-phicoef*F is the porosity. % % % % The boundary conditions are % Variable LS RS FE RE % Vf NO FLUX \nabla \cdot Vf = 0 % p FREE % Vf*G-D\nabla*G -F*(Vp) 0 % clear hmin clear fem aFw % Constants and domain data th0=pi/8; a=2.6; % semi axis of initial ellipse b=0.5; a1=1.0; b1=0.6; db=1.0; phicoeff=0.1; % phi=1-phicoeff*F; ht0=1; % cell height at the front edge hmin=0.120; % size of the mesh triangle xscale=1; % scale of dilation in x-direction arscale = 1.0; % arrow scale for the water velocity plot arrxnum=30; % number of arrows in x direction arrynum=25; % number of arrows in y direction % domain geometry kcheck = exist('hmin'); clear s c p %% %% Domain definition %% % upper ellipse ellipse1=ellip2(0,0,a,a,0); ellipse2=ellip2(0,-0.0*b,a,b,0); rec=rect2(-a,a,-a,0); ro=ellipse1-rec-ellipse2; % This removes all internal boundaries r=geomdel(ro); % r=ellip2(0,0,1.5,0.5,0); % Define objects objs={r}; names={'r'}; s.objs=objs; s.name=names; % % Curve objects % objs={}; names={}; c.objs=objs; c.name=names; objs={}; names={}; p.objs=objs; p.name=names; drawstruct=struct('s',s,'c',c,'p',p); fem.draw=drawstruct; fem.geom=geomcsg(fem); % This enables to see the numbering of the boundary % segments of the domain. %%%%% geomplot(fem, 'edgelabels', 'on', 'edgearrow', 'on'); pause; %%%%% % It appears that the following indices are given: % % % RE = 1 & 3 % FE = 2 & 4 % % Define aFw mode (F-actin + water) aFw.sdim={'x' 'y'}; aFw.dim = {'u','v','p'}; aFw.const = {'aa', a, 'bb', b, 'hh', ht0, ... 'V0x',0.0,'V0y',-0.25, ... 'Vpmin',0.1,'Vpmax',1.0,'Fmax',0.8,'Fmin',0.36, ... 'tt',th0,'phic', phicoeff}; aFw.expr = {'ksi0' 1.6 'ht' '1' ... '3*hh/2-hh/(2*bb)*y' ... 's0' 2.40 'F0' 1 ... 'Vpx' 'V0x-x*vel0(x,y,aa,bb,Vpmax,Vpmin,tt)/sqrt(x*x+y*y)' ... 'Vpy' 'V0y-y*vel0(x,y,aa,bb,Vpmax,Vpmin,tt)/sqrt(x*x+y*y)' ... ... 'F' 'Fmin+ggauss(x,y,0,aa,aa,0.6*aa,0.65*Fmax)' ... 'F' '(Fmax*(y-bb)+(aa-y)*Fmin)/(aa-bb)' ... }; aFw.geom=fem.geom; % Make initial mesh if kcheck aFw.mesh=meshinit(aFw,'hmax',hmin); else aFw.mesh=meshinit(aFw); end % aFw.bnd.ind = [3 1 2 2 1 3]; aFw.bnd.ind = [2 1 2 1]; aFw.shape = {shlag(1,'p') shlag(2,'u') shlag(2,'v')}; aFw.form = 'coefficient'; aFw.equ.a = {{'F*F' 0 0; ... 0 'F*F' 0; ... 0 ...'femval_stat(x,y,''-phic*Fx*ht+(1-phic*F)*htx'',''aF'')' ... 0 ...'femval_stat(x,y,''-phic*Fy*ht+(1-phic*F)*hty'',''aF'')' ... 0 ... }}; % F*Vf aFw.equ.be = {{ ... {0 0} {0 0} {'1/ksi0' 0}; ... {0 0} {0 0} {0 '1/ksi0'}; ... {0 0} {0 0} {0 0} ... }}; % 1/ksi0 \nabla p = \beta \cdot \nabla u aFw.equ.al = {{ ... {0 0} {0 0} {0 0}; ... {0 0} {0 0} {0 0}; ... {'-ht*(1-phic*F)' 0} {0 '-ht*(1-phic*F)'} {0 0} ... % {-1 0} {0 -1} {0 0} ... %simple case }}; % h*ph*[\nabla \cdot (Vf)] = -\nabla \cdot \alpha u aFw.equ.f={{... 'F*Vpx' ... 'F*Vpy' ... ... '-phic*0*F*(Vxx+Vyy)-phic*0*(Fx*Vpx+Fy*Vpy)' ... '-phic*0' ... }}; % No-flux BC on FE, RS and LS, and zero velocity on RE aFw.bnd.q = {... {'dnx' 'dny' '0'; ... 'dnx' 'dny' '0'; ... '0' '0' '0' }... {'0' '0' '0'; ... '0' '0' '0'; ... '0' '0' '0'} ... % {'0' '0' '0'; ... % '0' '0' '0'; ... % '0' '0' '0'} ... }; % aFw.bnd.g = {{'0' ... 'u*dnx+v*dny' ... % '0' ... 'u*dnx+v*dny' ... % '0'}}; aFw.bnd.h={... {'0' '0' '0'; ... '0' '0' '0'; ... '0' '0' '0'} ... {'1' '0' '0'; ... '0' '1' '0'; ... '0' '0' '0'}... % {'1' '0' '0'; ... % '0' '1' '0'; ... % '0' '0' '0'}... }; aFw.bnd.r={... {'0' '0' '0'} ... {'0' '0' '0'} ... % {'-Vout' '-Vout' '0'} ... }; aFw.xmesh = meshextend(aFw); % Solve the equations and plot the result. aFw.sol = femlin(aFw); postplot(aFw,'tridata','F','tribar','on',... 'Title','F-actin velocity',... 'arrowdata',{'Vpx' 'Vpy'},'arrowcolor','black'); %'arrowstyle','normalized'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); pause; % Plot the pressure and show the water velocity field on the same % picture. postplot(aFw,'tridata','p','tribar','on',... 'Title','pressure + water velocity',... 'arrowdata',{'u' 'v'},'arrowcolor','black','arrowscale',arscale,... 'arrowxspacing',arrxnum,'arrowyspacing',arrynum); %'arrowstyle','normalized'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); pause; postplot(aFw,'tridata','F','tribar','on',... 'Title','water flow + F-actin',... 'flowdata',{'u' 'v'},'flowlines',50,'flownormal','on','flowrefine',3); %'arrowstyle','normalized'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); pause; % Test the result by computing (\nabla*p-Vp)/(ksi*F) % postplot(aFw,'tridata','p','tribar','on',... % 'Title','pressure + water velocity check',... % 'arrowdata',... % {'femval_stat(x,y,''Vpx'',''aF'')-px/ksi0/femval_stat(x,y,''F'',''aF'')' ... % 'femval_stat(x,y,''Vpy'',''aF'')-py/ksi0/femval_stat(x,y,''F'',''aF'')'},... % 'arrowcolor','black','arrowscale',arscale); % daspect([1 xscale 1]); pbaspect([xscale 1 1]); % pause; % postcrossplot(aFw, 1, [1 2 5 6 3 4],'lindata','sqrt(v*v+u*u)','cont','on','lincolor','b',... % 'Title','velocity module'); % pause; % Define the fem mode for solution of the G-actin problem fem.sdim={'x' 'y'}; % Constants and domain data fem.expr = {'D' 4.0 'k' 5 'A0' 0.3 'B0' 0.1}; fem.dim = {'A','B'}; fem.bnd.ind = [2 1 2 1]; fem.shape = 2; % The coefficient form fem.form='coefficient'; % % The trick is not to use the original form of the equation with % the source term s0*F. Instead of it we replace this term by the % divergence of some flux term, namely s0*F=-\nabla \cdot (Vp*F). % % This trick leads to resolving of the positive source term problem % which otherwise would lead to the solution growing in time. fem.equ.f = {{'0' '0'}}; % {0, 0} in RHS fem.equ.c = {{'D' '0'; '0' 'D'}}; % -D \nabla^2 {a ,b} -> % \nabla \cdot (-c \nabla U) fem.equ.al = {{ ... {'femval_stat(x,y,''-u'',''aFw'')' 'femval_stat(x,y,''-v'',''aFw'')'} {0 0}; ... \% \nabla \cdot (Vf*{a, b}) {0 0} {'femval_stat(x,y,''-u'',''aFw'')' 'femval_stat(x,y,''-v'',''aFw'')'}... % \nabla \cdot (-alpha u) }}; fem.equ.ga = {{ ... {'femval_stat(x,y,''Vpx*F'',''aFw'')' 'femval_stat(x,y,''Vpy*F'',''aFw'')'}; ... {'0' '0'} ... }}; fem.equ.a={{'k' '-k'; '-k' 'k'}}; % {k(a-b), -k(a-b)} -> a U % % BC is \n \cdot (c \nabla U + alpha U - gamma) = g % % fem.bnd.q={... % '0'; ... % '0'; ... % '0' ... % }; % fem.bnd.g={ ... % {'0' '0'} ... % {'femval_statBnd(dom,s,''F*(Vpx*dnx+Vpy*dny)'',''aFw'')' '0'} ... % LS and RS % ... {'0' '0'} ... % RE % % }; fem.equ.init={{'A0' 'B0'}}; % fem.equ.da = 0.1; fem.mesh=meshinit(fem,'hmax',hmin); %fem=femdiff(fem); % For general form only fem.xmesh = meshextend(fem); % Solve the equations and plot the result. fem.sol = femlin(fem); % fem.sol=femtime(fem,'tlist',linspace(0,2.0,12),'report','on'); postplot(fem,'tridata','A','tribar','on',... 'Title','G-actin-profilin'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); axis off; pause; postplot(fem,'tridata','B','tribar','on',... 'Title','G-actin-thymosin'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); axis off; pause; % postplot(fem,'tridata',... % 'femval_stat(x,y,''ux+vy'',''aFw'')*A +(femval_stat(x,y,''u'',''aFw'')*Ax+femval_stat(x,y,''v'',''aFw'')*Ay)-D*(Axx+Ayy)-femval_stat(x,y,''s0*F'',''aF'')-k*(A-B)',... % 'tribar','on',... % 'Title','Error G-profilin'); axis equal; % pause; % postplot(fem,'tridata',... % 'femval_stat(x,y,''ux+vy'',''aFw'')*B +(femval_stat(x,y,''u'',''aFw'')*Bx+femval_stat(x,y,''v'',''aFw'')*By)-D*(Bxx+Byy)+k*(A-B)',... % 'tribar','on',... % 'Title','Error G-thymosin'); axis equal; % pause; postplot(fem,... 'tridata','A+B', ... ... 'A+B+femval_stat(x,y,''F'',''aFw'')',... 'tribar','on',... 'Title','Total G-actin'); daspect([1 xscale 1]); pbaspect([xscale 1 1]); % axis equal; pause; % plot of the profilin density on the front edge postcrossplot(fem, 1, [2 4], 'lindata','A','cont','on','lincolor','r'); hold on; % plot of the thymosin density on the front edge postcrossplot(fem, 1, [2 4],'lindata','B','cont','on','lincolor','b',... 'Title','G-actin'); hold off;