Updates for virus dynamics

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
+

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')
+

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));

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);
+

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
+
+
+