Matlab code for operatorial model of tumoral cell proliferation

hamiltoniansystem.m

function H=hamiltoniansystem(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu)
%Construct the full (time dependent) hamiltonian

NumTot=sum(Num);


H=sparse(NumTot,NumTot);
H=omega(1)*n{1}+omega(2)*n{2}+omega(3)*n{3}+...
    +mu(1)*(Hms)+...
        +mu(2)*(HmM)+...
          mu(3)*HsS+...
          mu(4)*(Hmg1)+...
          mu(5)*(Hmg2);


  

init.m

Num(1)=50; %number of bosonic modes healthy cells
Num(2)=150; %number of bosonic modes sick cells
Num(3)=2;   %number of bosonic modes med treat.
NumTot=sum(Num);

%Vacua vectors

for j=1:3
    phi0{j}=sparse(Num(j),1);
phi0{j}(1,1)=1;
end
phiG=sparse(NumTot,1);
phiG=sparse(kron(phi0{1},kron(phi0{2},phi0{3})));





%Pauli Matrixes
sp=[0 1;0 0];
sz=[1 0;0 -1];
Id=[1 0;0 1];
%Annihilation and Creation operators

for j=1:2
a{j}=sparse(Num(j),Num(j));
end

for h=1:Num(1)-1
    a{1}(h,h+1)=sqrt(h);
end
for h=1:Num(2)-1
    a{2}(h,h+1)=sqrt(h);
end

a{3}=sp;


p{1}=sparse(sparse(kron(a{1},sparse(kron(eye(Num(2)),eye(2))))));
p{2}=sparse(sparse(kron(eye(Num(1)),sparse(kron(a{2},eye(2))))));
p{3}=sparse(sparse(kron(eye(Num(1)),sparse(kron(eye(Num(2)),a{3})))));

%Number operators
for j=1:3
n{j}=sparse(p{j}'*p{j});
end

%Projector PS
j=0;
PSmax=sparse(Num(1)*Num(2)*Num(3),Num(1)*Num(2)*Num(3));
for j=0:Num(1)-1
    qq=(p{1}')^j*p{3}'*(p{2}'^(Num(2)-1))*phiG/(sqrt(factorial(j))*sqrt(factorial(Num(2)-1))); 
    qq2=(p{1}')^j*(p{2}'^(Num(2)-1))*phiG/(sqrt(factorial(j))*sqrt(factorial(Num(2)-1)));
    
PSmax=PSmax+outerproduct(qq,qq'.')+outerproduct(qq2,qq2'.');
end

%Main hamiltonian Operators
Hms=sparse(p{1}*p{2}'*n{1}+sqrt(Num(2))*p{1}*n{1}*PSmax);
HmM=sparse(p{2}'*n{2});
Hmg1=sparse(p{2}*n{2}*p{3});
Hmg2=sparse(p{2}*n{2}*n{3});
HsS=sparse(p{1}'*n{1});


%Inizialt state
phiin=(p{1}')^40*(p{2}')^10*p{3}'*phiG;
phiin=phiin/norm(phiin);

%Compute Evolution
tspan=[0:0.01:5];
[psi,T]=shroedingerevol(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu,phiin,tspan);


%Compute mean numbers%
for j=1:length(tspan)
nn(j)=psi(j,:)'.'*n{1}*psi(j,:).'/norm(psi(j,:))^2;
nn2(j)=psi(j,:)'.'*n{2}*psi(j,:).'/norm(psi(j,:))^2;   
nn3(j)=psi(j,:)'.'*n{3}*psi(j,:).'/norm(psi(j,:))^2;   
end

%mean values of hamiltonain operators%
h5=0;
for j=1:length(tspan)

mHmM(j)=psi(j,:)'.'*(mu(2)*exp(-((Num(2)-1))/(nn(2)))*HmM)*psi(j,:).'/norm(psi(j,:))^2;   
mHsS(j)=psi(j,:)'.'*(mu(3)*((Num(2)-1-real(nn2(j)))/(Num(2)-1))*HsS)*psi(j,:).'/norm(psi(j,:))^2; 
t=tspan(j);
mHmg1(j)=psi(j,:)'.'*(mu(4)*Hmg1)*psi(j,:).'/norm(psi(j,:))^2;   
mHms(j)=psi(j,:)'.'*(mu(1)*exp(-((Num(2)-1))/(nn(2)))*Hms)*psi(j,:).'/norm(psi(j,:))^2;
end


%Projectros of healthy state and sick state only%
for j=1:Num(1)
phis{j}=(p{1}')^j*p{3}'*phiG/(sqrt(factorial(j)))+(p{1}')^j*phiG/(sqrt(factorial(j)));
end
for j=1:Num(2)
phim{j}=(p{2}')^j*p{3}'*phiG/(sqrt(factorial(j)))+(p{2}')^j*phiG/(sqrt(factorial(j)));
end


for j=1:length(tspan)
    psin=psi(j,:)/norm(psi(j,:));
    probm(j)=0;
     probs(j)=0;
for k=1:Num(2)
probm(j)=probm(j)+abs(psin'.'*phim{k}).^2;
end
for k=1:Num(1)
probs(j)=probs(j)+abs(psin'.'*phis{k}).^2;
end
end

shroedingerevol.m

function [psi,T]=shroedingerevol(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu,phiin,tspan)

%Shroedinger evolution with RK45 code

options = odeset('RelTol',1e-4,'AbsTol',1e-4*ones(length(phiin),1))';
psi(length(tspan),length(phiin))=0;
T=tspan;
[T,psi] = ode45(@(t, z) secmem(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu,z,t),tspan,phiin,options);

 
 
 
 
 
 
    function jout=secmem(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu,z,t)
        
     
    t
        for j=1:3
        nn(j)=abs(z'*n{j}*z)/norm(z)^2;
        end
   
        mu(3)=mu(3)*((Num(2)-1-real(nn(2)))/(Num(2)-1))^1;
        mu(5)=mu(5);
        mu(2)=mu(2)*exp(-((Num(2)-1))/(nn(2)));
        mu(1)=mu(1)*exp(-((Num(2)-1))/(nn(2)));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
 
         
        H=hamiltoniansystem(Num,n,Hms,HmM,Hmg1,Hmg2,HsS,omega,mu);
        jout=-1i*H*z;