diff --git a/antigenN.m b/antigenN.m new file mode 100755 index 0000000..ad9124c --- /dev/null +++ b/antigenN.m @@ -0,0 +1,88 @@ +% function omegout = antigenN(node,gext) +% modified from coupleN.m on 2-22-12 + + +function [runv runx] = antigenN(N) + +if N == 1 + rtemp = 2; +else + rtemp = randpoisson(0.1, 1,N); +end +rtemp(N) = rtemp(N) + 2; +rg(1:N) = rtemp; + + +next = 1; +tnow = 0; +ypop(1) = 0; +for loop = 1:N + current_num = loop; + tstart = tnow ; + tend = tstart + 100*rg(loop); + tspan = [tstart tend]; + + ypop(2*loop) = 0.002; + ypop(2*loop+1) = 0; + + %ypop + + [t,y] = ode45(@f5,tspan,ypop); + [sy,nn] = size(y); + + nvirus = nn/2; + + for vloop = 1:nvirus + runv(next:(next+sy-1),vloop) = y(:,2*vloop); + runx(next:(next+sy-1),vloop) = y(:,2*vloop+1); + + ypop = y(sy,:); + + end + + next = next + sy; + tnow = tstart; +end + + + + function yd = f5(t,y) + + r = 2.4; % 2.4 + a = 2; % 2 + b = 0.1; % 0.1 0.1 0.01 + c = 1; % 0.1 1 1 + k = 1; % 0 1 1 % cross-reactive + q = 2.4; % 0 2.4 2.4 % cross-reactive + u = 0.2; % 0 0 1 % anti-immune + + vv = 0; + for nloop = 1:current_num + vv = vv + y(2*nloop); + end + + for nloop = 1:current_num + + v(nloop) = y(2*nloop); + x(nloop) = y(2*nloop+1); + z = y(1); + + dv(nloop) = v(nloop)*(r - a*x(nloop) - q*z); + dx(nloop) = -b*x(nloop) + c*v(nloop) - u *vv*x(nloop); + dz = k*vv - b*z - u*vv*z; + + end + + + + for omloop = 1:current_num + yd(2*omloop,1) = dv(omloop); + yd(2*omloop+1,1) = dx(omloop); + yd(1,1) = dz; + end + + end % end f5 + + +end % end antigenN + diff --git a/antigenNdriver.m b/antigenNdriver.m new file mode 100644 index 0000000..998f7aa --- /dev/null +++ b/antigenNdriver.m @@ -0,0 +1,40 @@ +% antigenNdriver.m +% 2-22-12 + + +clear +close all + + +N = 20; + +[runv runx] = antigenN(N); + +[sy,sx] = size(runv); + +xx = 1:sy; + +% figure(1) +% plot(x,runv(:,1),x,runx(:,1),x,runv(:,2),x,runx(:,2),x,runv(:,3),x,runx(:,3)) + +figure(1) +hold on +for loop = 1:N + plot(xx,runv(:,loop),'--',xx,runx(:,loop)) +end +legend('Virus','Immune') +hold off + + +v = threshmask(runv,0.001).*(runv); +x = threshmask(runx,0.001).*(runx); + + +smv = sum(v,2); +smx = sum(x,2); + + +figure(2) +plot(xx,smv,'b',xx,smx,'r') +legend('Virus','Immune') + diff --git a/fermi.m b/fermi.m new file mode 100755 index 0000000..abbd85b --- /dev/null +++ b/fermi.m @@ -0,0 +1 @@ +% Fermi function. % Usage: Fermi(x), where x is a vector. % y = 1./(1+exp(x)) function y = fermi(x) y = 1./(1+exp(x)); \ No newline at end of file diff --git a/randpoisson.m b/randpoisson.m new file mode 100644 index 0000000..91421d4 --- /dev/null +++ b/randpoisson.m @@ -0,0 +1,9 @@ +% function y = randpoisson(av,M,N) +% Creates M by N array with Poisson distribution of std = av + +function y = randpoisson(av,M,N) + +x = rand(M,N); + +y = -av*log(x); + diff --git a/threshmask.m b/threshmask.m new file mode 100644 index 0000000..1222955 --- /dev/null +++ b/threshmask.m @@ -0,0 +1,17 @@ +% function M = threshmask(A,thrsh) +% (Questionable functionality) + +function M = threshmask(A,thrsh) + +if thrsh == 0 + mxA = max(max(A)); + M = round(fermi(500*A/mxA)); +else +B = sign(thrsh).*1e6*(1 - A/(thrsh + 1e-20)); + +M = round(fermi(B)); + +end + + +