function [htimecourse,hcyclestats] = ccmodel(systemsize, ismutant)

if ~exist('ismutant')
    ismutant = false;
end

% Kinetic Parameters for wild type

[ k1,kp2,kpp2,kppp2, ...
  kp3,kpp3,J3,kp4,k4,J4, ...
  kp5,kpp5,k6,J5,k7,k8,J7,J8, ...
  k9,k10,J9,J10, ...
  k11,k12,kp12,kpp12,Kdiss, ...
  k13,k14, ...
  k15,kp16,kpp16,J15,J16, ...
  Vawee,Viwee,Jawee,Jiwee, ...
  Va25,Vi25,Ja25,Ji25, ...
  kpwee,kppwee,kp25,kpp25, ...
                        mu] = deal(...
  .03,.03,1,.1, ...
  1,10,.01,2,35,.01, ...
  .005,.3,.1,.3,1,.25,.001,.001, ...
  .1,.04,.01,.01, ...
  .1,.01,1,3,.001, ...
  .1,.1, ...
  1.5,1,2,.01,.01, ...
  .25,1,.01,.01, ...
  1,.25,.01,.01, ...
  .15,1.3,.05,5, ...
                    0.005);
                
%  Kinetic parameters for double-mutant type

if ismutant
    [kppwee,kp25,kpp25] = deal(0.3,.02,.02);
end

% Initial consitions

% [Cdc13T, preMPF, Ste9, Slp1T, Slp1, IEP, Rum1T, SK, M]

init = [1.5, .05, 0, 0, 0, 0, .01, 0, 1.8]';

% init = [0.6, 0.018, 0.6, 2.13, 2.13, 0.925, 0.04, 0.535, 1.99]';

% init(9,1) = 1.99;
% 
% init(1) = 0.59999999999999998 ; 
% init(2) = 0.018045828316532363 ; 
% init(3) = 0.59999999999999998 ;
% init(4) = 2.1299999999999999 ; 
% init(5) = 2.1299999999999999 ; 
% init(6) = 0.92500000000000004 ; 
% init(7) = 0.040000000000000001 ;
% init(8) = 0.53500000000000003 ; 



omegainv = 1./systemsize;

zeros8x8 = zeros(8);
ones8x8 = ones(8);

ind = find(tril(ones8x8));
indl = find(tril(ones8x8,-1));
indr = find(triu(ones8x8,1));
[indlx,indly] = ind2sub([8 8], indl);
indltrans = sub2ind([8 8], indly, indlx);
indmat = zeros8x8;
indmat(ind) = 10:45;
indiag = diag(indmat);

% compind = repmat([1:8]',1,2);
% S = repmat([1 -1],8,1);

% Goldbeter-Koshland function

% G = @(a,b,c,d)
% 2*a.*d./(b-a+b.*c+a.*d+sqrt((b-a+b.*c+a.*d).^2-4*a.*d.*(b-a)));

htimecourse = @timecourse;
hcyclestats = @cyclestats;


% solve the differential equations  


function [ct,bm,dm] = cyclestats(ncells,ncycles)
    
    tfinal = 6e10;
        
%         mu_bm = 9.9/8.2;
%         sd_bm = 1.96/8.2;        
%         init_mass = normrnd(mu_bm,sd_bm,ncells,1);
        
%         init_mass = [11.3 10.4 9.4 5.7]./8.2;

        if ncells==1
            init_mass = init(9);
        else
            init_mass = 1.8*randsample(linspace(0.48,2.4,1e3).',ncells);
        end
        
        indnonnegative = [1:9 indiag'];
        y0 = [init;zeros(36,1)];

        odefun = @odemes;
        
        options = odeset( 'Events', @events, ...
            'RelTol', 1e-4, 'InitialStep',0.0005, 'MaxStep', 20, ...
            'NonNegative', indnonnegative, 'MaxOrder', 2, 'BDF', 'off' );

        
        ncycles = ncycles + 1;
        tete = zeros(ncells,ncycles);
        meme = zeros(ncells,ncycles);
        
        for i=1:ncells
            tstart = 0;
            y0(9) = init_mass(i);
            for j=1:ncycles
                [t,y,te,ye] = ode15s(odefun,[tstart tfinal],y0,options);
                nt = length(t);
                tete(i,j) = te;
                meme(i,j) = ye(9);
                y0 = y(nt,:).';
                y0(9) = .5*y0(9);
                tstart = t(nt);
                clear t y te ye                
            end
        end
        
        ct = diff(tete, 1, 2);
        ct = ct(:);
        dm = meme(:,2:end);
        dm = dm(:);
        bm = meme(:,1:end-1);
        bm = 0.5*bm(:);       
end 

function solout = timecourse(tfinal,ismesoscopic) 
        
        if ismesoscopic
            indnonnegative = [1:9 indiag'];
            y0 = [init;zeros(36,1)];
            odefun = @odemes;
        else
            indnonnegative = 1:9;
            y0 = init;
            odefun = @odemac;
        end

        options = odeset( 'Events', @events, ...
            'RelTol', 1e-4, 'InitialStep',0.0005, 'MaxStep', 20, ...
            'NonNegative', indnonnegative, 'MaxOrder', 2, 'BDF', 'off' );
        refine = 1;
        
        tstart = 0;

        tt = tstart;
        yy = y0.';


        while tstart<tfinal
            % Solve until the first terminal event.
            [t,y] = ode15s(odefun,[tstart tfinal],y0,options);
            % Accumulate output.
            nt = length(t);
            tt = [tt; t(2:nt)];
            yy = [yy; y(2:nt,:)];
            
            y0 = y(nt,:).';
            y0(9) = .5*y0(9);
            tstart = t(nt);
            clear t y
        end
        
        solout.t = tt;                
        M = yy(:,9).';
        u = yy(:,1:8).';
        [u(9,:),ktr] = getMPFTrimer(u);
        ktf = G(k15*M,kp16+kpp16*u(9,:),J15,J16);
        [u,M,ktr,ktf] = deal(u',M',ktr',ktf');
        solout.y = [u, M, ktr, ktf];
        
        if ismesoscopic
            k = zeros8x8;
            k(ind) = 1:36;            
%             indoutl9 = find(tril(ones(9),-1));
            indout = find(tril(ones(9)));
            kout = zeros(9);
            kout(indout) = 1:45;
%             indiag9 = diag(kout);
            w = yy(:,10:end);
            m(:,9) = (u(:,1)-ktr).^2 + w(:,k(1,1)) - u(:,9).^2 - ...
                2*u(:,1) + 2*(2-ktr./u(:,1)).*ktr.*u(:,2) + ...
                2*((ktr.^2+6*ktr.*u(:,2))./u(:,1).^2 - ...
                4*ktr.^2.*u(:,2)./u(:,1).^3 - 1).*w(:,k(2,1)) + ...
                (1-ktr./u(:,1)).^2.*(u(:,2).^2+w(:,k(2,2))) + ...
                ((3*ktr./u(:,1)).^2-2).*ktr.*u(:,2).^2.*w(:,k(1,1))./u(:,1).^3;

            m(m(:,9)<0,9) = 0;
            m(:,1) = u(:,1).^2 + w(:,k(1,1)) - w(:,k(2,1)) - ...
                (u(:,2)+u(:,9)).*u(:,1) - ktr.*(u(:,1)-u(:,2));
            
            for i=2:8
                m(:,i) = w(:,k(i,1)) + ...
                 (u(:,1)-u(:,2)-u(:,9)-ktr+3*ktr.*u(:,2)./u(:,1)).*u(:,i)+ ...
                                       (3*ktr./u(:,1)-1).*w(:,k(i,2));
            end

            wm(:,kout(end,:)) = m;
            wm(:,nonzeros(tril(kout(1:end-1,1:end-1)))) = w;

            solout.w = wm;
            
%             [indx,indy] = ind2sub([9 9], indoutl9);
            [indx,indy] = ind2sub([9 9], indout);
            wm = wm./(u(:,indx).*u(:,indy));            
            solout.nsr = wm;
            solout.cv = sqrt(wm(:,diag(kout)));
            solout.kout = kout;
%             solout.cv = sqrt(var9)./u;            
            
            
%             var9 = wm(:,indiag9);
%             solout.cv = sqrt(var9)./u;
%             size(wm)
%             max(indl9(:))
%             size(var9)            
%             max(indx(:))
%             max(indy(:))
            
%             indl9 = nonzeros(tril(kout,-1));
%             wm(:,indl9) = wm(:,indl9)./sqrt(var9(:,indx).*var9(:,indy));
%             wm(:,indiag9) = 0;
%             solout.cc = wm;        
            
        end
%         
    end 

 % Macroscopic ODE model

    function dudt = odemac(t,u)

        MPF = getMPFTrimer(u);
        ktf = G(k15*u(9),kp16+kpp16*MPF,J15,J16);
        kwee = kpwee + (kppwee-kpwee)*G(Vawee,Viwee*MPF,Jawee,Jiwee);
        k25 = kp25 + (kpp25-kp25)*G(Va25*MPF,Vi25,Ja25,Ji25);

        kpu = kp2 + kpp2*u(3) + kppp2*u(5);

        JInv = 1./[J4+u(3); J3+1-u(3); J7+u(4)-u(5); ...
                                        J8+u(5); J9+1-u(6); J10+u(6)];

        v(8,:) = [k13*ktf, k14*u(8)];
        v(1,:) = [k1*u(9), kpu.*u(1)];
        v(2,:) = [kwee*(u(1)-u(2)), (k25+kpu).*u(2)];
        v(3,1) = (kp3+kpp3*u(5)).*(1-u(3)).*JInv(2);
        v(3,2) = (kp4*u(8) + k4*MPF).*u(3).*JInv(1);
        v(4,:) = [kp5 + kpp5 * MPF^4/(J4^4+MPF^4), k6*u(4)];
        v(5,1) = k7*u(6).*(u(4)-u(5)).*JInv(3);
        v(5,2) = (k6+k8*JInv(4))*u(5);
        v(6,:) = [k9*MPF*(1-u(6)).*JInv(5), k10*u(6).*JInv(6)];
        v(7,:) = [k11, (k12+kp12*u(8)+kpp12*MPF).*u(7)];

        %----------------------------------------------------------------
        
        dudt = [v(:,1) - v(:,2); mu*u(9)];

    end

% Mesoscopic ODE model

    function dydt = odemes(t,y)
        
        % Quantities in pseudo-steady-state
        
        MPF = getMPFTrimer(y);
        ktf = G(k15*y(9),kp16+kpp16*MPF,J15,J16);
        kwee = kpwee + (kppwee-kpwee)*G(Vawee,Viwee*MPF,Jawee,Jiwee);
        k25 = kp25 + (kpp25-kp25)*G(Va25*MPF,Vi25,Ja25,Ji25);
        
        % Tempory variables
        
        kpy = kp2 + kpp2*y(3) + kppp2*y(5);
        kpMPF = kp4*y(8) + k4*MPF;
        Jinv(:,3) = 1./[J4+y(3); J3+1-y(3); J7+y(4)-y(5); ...
                                   J8+y(5); J9+1-y(6); J10+y(6)];
        Jinv(:,1) = Jinv(:,3);
        Jinv(:,2) = Jinv(:,1).^2;
        Jinv(:,3) = Jinv(:,1).*Jinv(:,2);

        % Coefficients of Flux and Diffusion
        
        v(8,:) = [k13*ktf, k14*y(8)];
        v(1,:) = [k1*y(9), kpy.*y(1)];
        v(2,:) = [kwee*(y(1)-y(2)), (k25+kpy).*y(2)];
        v(3,1) = (kp3+kpp3*y(5)).*(1-y(3)).*Jinv(2,1);
        v(3,2) = kpMPF.*y(3).*Jinv(1,1);
        v(4,:) = [kp5 + kpp5 * (MPF.^4)./(J4^4+MPF.^4), k6*y(4)];
        v(5,1) = k7*y(6).*(y(4)-y(5)).*Jinv(3,1);
        v(5,2) = (k6+k8*Jinv(4,1)).*y(5);
        v(6,:) = [k9*MPF*(1-y(6)).*Jinv(5,1), k10*y(6).*Jinv(6,1)];
        v(7,:) = [k11, (k12+kp12*y(8)+kpp12*MPF).*y(7)];
        
        
        f = v(:,1) - v(:,2);
        b = v(:,1) + v(:,2);
        
        % Jacobian
        
        A = zeros8x8;
        A(1,[1 3 5]) = [-kp2-kpp2*y(3)-kppp2*y(5), -kpp2*y(1), -kppp2*y(1)];
        A(2,[1:3 5]) = [kwee, -kwee-k25+A(1,1), -kpp2*y(2), -kppp2*y(2)];
        A(3,3) = -J4*(kpMPF).*Jinv(1,2) - J3*(kp3+kpp3*y(5)).*Jinv(2,2);
        A(3,[5 8]) = [kpp3*(1-y(3)).*Jinv(2,1), -kp4*y(3).*Jinv(1,1)];
        A(4,4) = -k6;
        A(5,4) = k7*J7*y(6).*Jinv(3,2);
        A(5,[5 6]) = [-k6-A(5,4)-k8*J8*Jinv(4,2), k7*(y(4)-y(5)).*Jinv(3,1)];
        A(6,6) = -k9*MPF*J9*Jinv(5,2) - k10*J10*Jinv(6,2);
        A(7,[7 8]) = [-k12+kpp2*MPF-kp12*y(8), -kp12*y(7)];
        A(8,8) = -k14;
        
        % Errors arising from Flux and Diffusion
        w = zeros8x8;
        w(ind) = y(10:end);
        w(indltrans) = w(indl);
        
        e(8,2) = 0;
        e(1,2) = kpp2*w(1,3) + kppp2*w(1,5);
        e(2,2) = kpp2*w(2,3) + kppp2*w(2,5);
        e(3,1) = -w(3,3)*J3*(kp3+kpp3*y(5)).*Jinv(2,3) ...
                                        - w(3,5)*kpp3*J3*Jinv(2,2);
        e(3,2) = -w(3,3)*J4*kpMPF.*Jinv(1,3) + w(3,8)*kp4*J4*Jinv(1,2);
        e(5,1) = (2*w(4,5)-w(4,4)-w(5,5))*k7*J7*y(6).*Jinv(3,3) ...
                                       + (w(4,6)-w(5,6))*k7*J7*Jinv(3,2);
        e(5,2) = -w(5,5)*k8*J8*Jinv(4,3);
        e(6,:) = -w(6,6)*[k9*MPF*J9*Jinv(5,3), k10*J10*Jinv(6,3)];
        e(7,2) = kp12*w(7,8) ;        
        
        ef = e(:,1) - e(:,2);
        eb = e(:,1) + e(:,2);        
                
        % Construct ODEs
        
        Aw = A*w;
        dwdt = tril(Aw + Aw.');
        dydt(10:45,1) = dwdt(ind);
        dydt(indiag) = dydt(indiag) + omegainv.*(b+eb);
        dydt(1:9) = [f + ef; mu*y(9)];
    end

% Check point for completion of one cell division

    function [value,isterminal,direction] = events(t,y)
        value = getMPFTrimer(y) - 0.1;
        isterminal = 1; % reset at the end of M-phase
        direction = -1;        
    end

% Compute MPF and ktr (assumed in pseudo-steady-state)

    function [MPF,ktr] = getMPFTrimer(y)
        Sigma = y(1,:) + y(7,:) + Kdiss;
        ktr = y(1,:).*y(7,:);
        ktr = 2*ktr./(Sigma+sqrt(Sigma.^2-4*ktr));
        if size(y,1)>9
            MPF = y(1,:) - y(2,:) - ktr + (ktr./y(1,:)).* ...
                        (y(2,:).*(1+y(10,:)./y(1,:).^2)-y(11,:)./y(1,:)) ;
        else
            MPF = (y(1,:)-y(2,:)).*(y(1,:)-ktr)./y(1,:);
        end
    end

% Goldbeter-Koshland function

    function gout = G(a,b,c,d)
        gout = 2*a.*d./(b-a+b.*c+a.*d+sqrt((b-a+b.*c+a.*d).^2-4*a.*d.*(b-a)));
    end

end