giordano2020

[anandijain]

@paulflang

dont merge

this script requires SciML/SBMLToolkit.jl#108

[anandijain]

we currently can recreate the plot in fig 2b. the weird thing is they scale everything down by 10e3 so when we plot infected we get it going to ~5e-4, which somehow is 60% of the population

here is the comparisons
from paper:

us:

[anandijain]

xcur=[200 20 1 2 0]'/60000000;
x=xcur;
L=[-r1(1)+alpha(1)+x(2)/x(1)*beta(1)+x(3)/x(1)*gamma(1)+x(4)/x(1)*delta(1);
    x(1)/x(2)*epsilon(1)-r2(1);
    x(1)/x(3)*zeta(1)-r3(1);
    x(2)/x(4)*eta(1)+x(3)/x(4)*theta(1)-r4(1); 0];
    
SUScur=(60000000-200-20-1-2)/60000000;
GU=0;
MO=7/60000000;
GUD=0;
GUDvac=0;
currentcases=0;

%%%%%%%%%%%% simulation 
for k=1:N_measures-1
        for i=(T_measures(k)-1)/T_camp+1:(T_measures(k+1)-1)/T_camp
               vacc(i)=vac(k); 
               
               if trigger==1
                  vacc(i)=max([(1-currentcases*60)*VACC(i,vvv) 0]);
                   %vacc(i)=max([(phi-currentcases*60*phi)*vac(k)/phi 0]);
                  end
               if giuseppe==1
                   vacc(i)=VACC(i,vvv);
               end
            if phiaffine==0
            SUS=SUScur-T_camp*SUScur*(CS(k,:)*xcur+vacc(i)); %\phi*S
            else
          SUS=SUScur-T_camp*SUScur*(CS(k,:)*xcur)-T_camp*vacc(i);% \phi
            end
          x=xcur+T_camp*(FF(:,:,k)+b*SUScur*CS(k,:))*xcur;
          L=inv(diag(xcur))*(FF(:,:,k)+b*SUScur*CS(k,:))*xcur;
          I(i)=x(1); 
          D(i)=x(2);
          A(i)=x(3);
          R(i)=x(4);
          T(i)=x(5);
          S(i)=SUS;
          NP(i)=epsilon(k)*I(i)+(theta(k)+mu(k))*A(i);
          R0cur(i)=R0(k);
          q(i)=alpha(k)+beta(k)*D(i)/I(i)+gamma(k)*A(i)/I(i)+delta(k)*R(i)/I(i);
          r1cur(i)=r1(k);
          YS(i)=CS(k,:)*x;
          if phiaffine==0
         
           YH(i)=CH(k,:)*x+vacc(i)*SUScur; %\phi*S
           YHDvac(i)=CHdiag(k,:)*x+vacc(i)*SUScur;
         
          else
              YH(i)=CH(k,:)*x+vacc(i);  %phi
               YHDvac(i)=CHdiag(k,:)*x+vacc(i);
          end
          YHD(i)=CHdiag(k,:)*x;
          YE(i)=CE(k,:)*x;
          GU=GU+T_camp*YH(i);
          GUD=GUD+T_camp*YHD(i);
          GUDvac=GUDvac+T_camp*YHDvac(i);
          MO=MO+T_camp*YE(i);
          H(i)=GU;
          HD(i)=GUD;
          HDvac(i)=GUDvac;
          E(i)=MO;
          SUScur=SUS;
          xcur=x;
          currentcases=R(i)+T(i)+D(i);
          LI(i)=L(1); 
          LD(i)=L(2);
          LA(i)=L(3);
          LR(i)=L(4);
          LT(i)=L(5);
          REFF(i)=(alpha(k)+beta(k)*x(2)/x(1)+gamma(k)*x(3)/x(1)+delta(k)*x(4)/x(1))*SUS/r1(k);
          REFFcur(i)=(alpha(k)+beta(k)*x(2)/x(1)+gamma(k)*x(3)/x(1)+delta(k)*x(4)/x(1))/r1(k);
          
    end
end

this is from simulation.m in the V model. if you squint really hard you can sort of see a model

[anandijain]

okay actually looking at the paper i can make the model, but not any of the data stuff. matlab code pretty unusable
V

first model

i dont know if this is in the paper, but its worth noting that they not only add V state, but an extra edge from R to E

should have this done in like 10

[anandijain]

default simulation of V model not looking too unreasonable. havent compared it with anything. but it looks like an SIR simulation 😅

[anandijain]

not sure what i did, but this model does show us peaking at 50 days with ~50 % of the populatiuon.

this is a lot more sane than the first model

[anandijain]

@ChrisRackauckas ready