
% 13/04/2022
% Nikhil Walani, UPC
% Main_File to simulate the combined effects of Protein diffusion, actin
% cortex drag and Branched actin push on budded membrane


clear all
close all
addpath Functionals:Interpolation_Geometry:Reparametrization:Input_Files:Matrices:Output_Files

global phi C_1 dt ini_Val Material Constraints V C_psi

% Preload the budded shape
input_file ='Input_Files/Budded_Membrane_Eq_F0_0_12.mat';
load(input_file)

% Options for the matlab solver
MaxFunEvals=3000;
MaxIter=50000;
tol_X=1e-3;

options = optimoptions('fmincon','Algorithm','interior-point' ,'Display','final', 'CheckGradient', false, ...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,'MaxFunctionEvaluations',MaxFunEvals);

options = optimoptions(options,'TolX',1e-7,'Display','final');

% Specify whether solving for an open patch of membrane 
iopen=ini_Val.iopen;

% Do you want to impose the fixed volume with support? 0: No 1: Yes
Constraints.iVol_s = 1;


%How many steps prior to final shape you wish to start
S_back =1;
x = x_hist(:,end-S_back+1);


%Constraints.C_s=C_s_eq(end-S_back);
dt = [];
time = [];

%% Parameters
%Membrane viscosity
Material.xi_L2 = 0;
%Membrane viscosity
Material.mu = 3.65e-3;
%Cortex viscosity (of a thin layer)
Material.eta_c = 5e3 / 1000 / 10  ;
%Membrane support friction
Material.eta_s = 2.4e-1 * 5     * 1;
%Parameters for z-cutoff of friction.
Material.tmp_s = 10;
Material.h_cut_off = 1;
% Force at edge of membrane
Constraints.f0 = -0.12 ;


% EFFECTIVE VISCOSITY/ELASTICITY OF THE REST OF THE SYSTEM
Material.eta_END = Material.eta_s * pi * x(end-1)^2 * log(100);
Material.xi_Cor_perp = 0;
Material.K_outer = 1 * 0;

% This is membrane I-BAR protein
ap = 3;
Material.Chi=0;
Material.RT_a0 = 1/20 / ap;
Material.phi_max=0.35;
Constraints.C0=-0.03  ;
Diff = 1; % this is 0.1 um^2/s
Material.eta_a0 =  1/20 / Diff / ap;
Material.phi_bar = 0.05;
phi_bar = Material.phi_bar;



% This is regulator
Material.D = 1 / 5000 ;
Material.k1= 1 / 5000;
Material.k2= 1 / 5000;
Material.enrichment_treshold = 2;

% Coefficients that transform psi into contractility as f1 - f2*psi
Material.f1 = 2;
Material.f2 = 0.7;


%% Times
t_f = 1000 * 30 * 4;
time = [0];
x_hist = x;
dt = 0.001 * 100 ;
i_max = 3000;
dt_max = 10 * 10;


%% Normal code

ini_Val.UU =ini_Val.U(4:end-3);

C_phi = phi_bar * ones(size(C_phi));
[phi, d_phi] = interpolate_func_2(Invariant,C_phi);
[phi_s, d_phi_s] = interpolate_func_2(Invariant_sample,C_phi);

C_psi = 0 * ones(size(C_phi));
[psi, d_psi] = interpolate_func_2(Invariant,C_psi);


[C,C1,C2] = interpolate(Invariant,x);
Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2,phi);


[C_s,C1_s,C2_s] = interpolate(Invariant_sample,x);
% Update the local area constraints
Constraints.S0_L = localArea_functional(Aux);
Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
[tS] = Area_functional(Aux);
[tV] = Volume_functional_support(Aux);

figure(2)
hold on
plot(C_s( :,1),C_s( :,2))
axis equal
hold on
plot(x(1:2:end),x(2:2:end),'ko-')

C_phi_hist = [C_phi];
phi_s_hist = [phi_s];
C_psi_hist = [C_psi];
C_s_hist = [C_s];
tot_ener_hist = [];
area_hist = [tS];
Chm_p_hist= [];
force_hist=[Constraints.f0];
force_hist_z=[];
Change_V0 = [];
sigma = [];

x_0=x;
x_1=x;
ct_rep=0;
dt_inc=0;
ini_Val.cc = 0.2;
ini_Val.ce = 1;
i=0;
epsl=1e-4;
ctr=0;
count=0;

Cortex_Flat_Updated(psi);
V_hist = [V];

% Force free velocity of growth of branched network
v0_dendrite = 9e-3; % Set this to 0 to consider drag by cortical actin
% Stall Force
F0          = 0.22;
% Governs the decay in growth velocity of branched network
decay_f     = 20;
% Specifies whether the dendritic network is in contact with the membrane
contact     = 0;
% Change in magnitude of height (just to initialize here)
Ht_Chng     = 0;

%Cortex velocity on cortex mesh
[v_c_c]=interpolate_func_2(Invariant,V);
%Cortex velocity on membrane
Invariant_mem_to_cortex=precal_generic(C(:,1)/Material.R0);
[v_c_m]=interpolate_func_2(Invariant_mem_to_cortex,V);


Proj_area_0 = pi * x(end-1)^2;
Areal_Strain = 1;
max_psi_hist = [];
max_phi_hist = [];
max_V_c_hist = [];
tip_hist       = [];
tip_force_hist = [];

Material.R_ini = x(end-1);

% Counter for number of time steps (inluding non-converged time stamps)
k = 0;

if(Material.phi_bar==0)
    t_c = 0;
else
    t_c = t_f; % time to start the counter for Arp2/3 actin
end
%%
while time(end) < t_f & i < i_max & k < 100000
    i = i+1;
    k = k+1;
    fprintf('Instant %10.5f towards %10.5f\n', time(end)+dt, t_f)

    % To set the start time for growth of branched network
    if(t_c==t_f && max(psi)>1e-2)
        t_c = time(end);
    end

    if(t_c < t_f && max(psi)<5e-3 && length(beq)==4)
        v0_dendrite = 0;
        beq(4)=[];
        Aeq(4,:)=[];
        contact = 0;
    end

    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);

    % To reparameterize
    max(abs(Aux.dau));
    if (any(abs(Aux.dau) > max_dau) || mod(i,20)==0)
        disp('**********************')
        disp('We reparametrize')
        [C_fit,phi_fit] = original_shape_to_fit(x, C_phi); % at points uniformly spaced in arc-length
        %[C_fit] = original_shape_to_fit(x); % at points uniformly spaced in arc-length
        ce=ini_Val.ce; cc=ini_Val.cc;
        %creparametrize(x,ce,cc);
        x_new = reparametrize(C_fit,x);
        x=x_new;
        % Update protein coverage
        C_phi_new = reparametrize_func(phi_fit,C_phi);
        C_phi=C_phi_new;
        [phi, d_phi, dd_phi] = interpolate_func_2(Invariant,C_phi);


        % Update the local area constraints
        [C,C1,C2] = interpolate(Invariant,x);
        Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2);
        Constraints.S0_L = localArea_functional(Aux);
        Constraints.S0_L = Constraints.S0_L * Constraints.S0 / sum(Constraints.S0_L);
        ct_rep=ct_rep+1;
        if(ct_rep==2)
            max_dau=max_dau*1.15;
            if(max_dau>200)
                max_dau = 200;
            end
            ct_rep=0;
        end
        pause
    else
        ct_rep=0;
    end

    % Solving for protein concentration given membrane shape
    phi_1 = phi;
    C_phi_1 = C_phi;
    if i>0
        C_phi_k = C_phi;
        [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
        R = Residual_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k);
        norm(R);
        M=matrix_Mass(Aux,C,C1,C2);
        A=matrix_AP2(Aux,C,C1,C2);

        ctr=1;
        while (norm(R) > 1e-10 && ctr <10)
            K=matrix_K_mod_Ener_Aug(Aux,C,C1,C2,phi_k);
            B=matrix_B_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k);
            JD=matrix_JD(Aux,C,C1,C2,phi_k,d_phi_k,dd_phi_k);
            J = (M + A + K + B +JD);
            J(:,end) = 0;
            J(end,:) = 0;
            J(end,end)=1;
            R(end)=0;

            C_phi_new = C_phi_k - J'\R;
            C_phi_k= C_phi_new;
            [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
            R = Residual_mod_Ener_Aug(Aux,C,C1,C2,phi_k,d_phi_k,phi_1,dd_phi_k);
            ctr=ctr+1;
        end
        [phi_k, d_phi_k, dd_phi_k] = interpolate_func_2(Invariant,C_phi_k);
        C_phi = C_phi_k;
        phi = phi_k;
        d_phi = d_phi_k;
        dd_phi = dd_phi_k;
    end


    % Map phi to the flat cortex
    [N ,EDGES ,BIN ] = histcounts(C(:,1),ini_Val.UU*Material.R0);
    [N_,EDGES_,BIN_] = histcounts(Invariant.gaussp,ini_Val.UU);
    phi_flat_ = zeros(size(Invariant.gaussp));
    for kk = 1:length(ini_Val.UU)-1
        if sum(BIN==kk)>0
            phi_flat(kk) = max(phi(BIN==kk));
        else
            phi_flat(kk)=phi_flat(kk-1);
        end
        phi_flat_(BIN_==kk) = phi_flat(kk);
    end

    % Solving for regulatory specie psi
    psi_1 = psi;
    C_psi_1 = C_psi;

    M_psi=matrix_Mass_psi;
    K_psi=matrix_Stiffness_psi;
    Rhs_psi = vector_rhs_psi(psi_1,phi);
    C_psi = (M_psi + K_psi)\Rhs_psi;
    [psi, d_psi] = interpolate_func_2(Invariant,C_psi);
    
    % Velocity of actin cortex depending on updated psi
    Cortex_Flat_Updated(psi);


    % Solving for membrane shape
    x_1 = x;
    [C_1,C1_1,C2_1] = interpolate(Invariant,x_1);
    [x, fval, exitflag, output, lambda] = fmincon(@Incremental_Mech_Rayleighian_S,x,[],[],Aeq,beq,[],[],@Nonlin_constr,options);



    % Adaptive time stepping to reduce the time step
    % for large changes in shape and increase for small changes
    if i>0
        if  max(abs(x-x_1)) > 0.5 || count > 5
            disp('*****************')
            disp('We reduce the timestep')
            dt = dt/2;
            if(contact)
                Ht_Chng  = Ht_Chng/2;
                beq(4)   = beq(4) - Ht_Chng;
            end
            pause(0.2)
            x = x_1;
            continue
        end
        if max(abs(x-x_1)) < 5e-2
            dt = dt * 1.5;
            if dt >= dt_max
                dt = dt_max;

            else
                disp('*****************')
                disp('We increase the timestep')
            end
        end

    end

    Material.eta_END = Material.eta_s * pi * x(end-1)^2 * log(100);

    [C,C1,C2] = interpolate(Invariant,x);
    Aux = auxiliary_calculations(Invariant,ini_Val.p,Constraints.C0,C,C1,C2,phi);

    Invariant_mem_to_cortex=precal_generic(C(:,1)/Material.R0);
    
    % Condition for establishment of contact between membrane and branched actin described
    % later. If contact established, then velocity of membrane apex is
    % determined by the force acting at the tip
    if(contact)
        disp(lambda.eqlin(4))
        Ht_Chng = dt*v0_dendrite*(exp(decay_f*lambda.eqlin(4)) - exp(-decay_f*F0));
        beq(4)  = beq(4) + Ht_Chng;
    end

    % Storage of data for each time stamp
    time = [time time(end)+dt];
    tip_hist = [tip_hist x(2)];
    
    if(contact)
        tip_force_hist = [tip_force_hist lambda.eqlin(4)];
    else
        tip_force_hist = [tip_force_hist 0];
    end

    C_psi_hist = [C_psi_hist C_psi];
    [phi_s] = interpolate_func_2(Invariant_sample,C_phi);
    phi_s_hist = [phi_s_hist phi_s];
    [C_s,C1_s,C2_s] = interpolate(Invariant_sample,x);
    C_s_hist = [C_s_hist C_s];
    V_hist = [V_hist V];
    area_hist = [area_hist Area_functional(Aux)];
    max_psi_hist  = [max_psi_hist  max(psi)];
    max_phi_hist  = [max_phi_hist  max(phi_s)/Material.phi_bar];


    % Obtain membrane tension field as Lagrange multiplier to local area
    % constraint
    if (Constraints.iSurf_loc ==1) %&& count==0
        tension = interpolate_func(Invariant,lambda.eqnonlin(2:end),Invariant.lower);


    else
        tension = lambda.eqnonlin*ones(size(Invariant.gaussp));
        %count=1;
    end
    sigma = [sigma tension];

    %plot with respect to arclength
    D= (Aux.au.* Invariant.gweight);
    S=zeros(size(D));
    for i=1:size(D)
        S(i)= sum(D(1:i));
    end
    %%%%%%%%%%%%%%%%%%%%%
    figure(3)
    hold off
    plot(C_s( :,1),C_s( :,2))
    hold on
    plot(x(1:2:end),x(2:2:end),'ko')
    axis equal
    
    %Cortex velocity on membrane
    [v_c_m]=interpolate_func_2(Invariant_mem_to_cortex,V);
    [v_c_c]=interpolate_func_2(Invariant,V);
    quiver(Invariant.gaussp*Material.R0,Invariant.gaussp*Material.R0*0,v_c_c',v_c_c'*0)
    tmp_z  = 0.5*(1-tanh(Material.tmp_s*(C(:,2)-Material.h_cut_off)));
    v_c_m = v_c_m .* tmp_z;
    quiver(C(:,1),C(:,2),v_c_m,v_c_m*0)

    max_V_c_hist  = [max_V_c_hist  max(v_c_c)];

    % Condition for establishment of contact between branched actin and
    % membrane 
    if(v0_dendrite*(1-exp(-decay_f*F0))*(time(end)-t_c)>x(2) && contact==0)
        dt          = dt/5;
        Aeq         = [Aeq; zeros(1,length(Aeq))];
        Aeq(4,2)    = 1;
        Ht_Chng     = dt*v0_dendrite*(1-exp(-decay_f*F0));
        beq(4)      = x(2) + Ht_Chng;
        contact     = 1;
        disp('We establish contact between growing dendritic network and membrane')
    end


    Proj_area_i = pi*x(end-1)^2;
    Areal_Strain = [Areal_Strain (pi*Material.R0^2 - Proj_area_i )/(pi*Material.R0^2 - Proj_area_0 )];
    figure (7)
    title('Normalized excess area vs Time')
    plot(time,Areal_Strain,'--ro','Markersize',4,'LineWidth',2)
    xlabel('Time','FontSize',20,'FontWeight','bold','Color','k')
    ylabel('Normalized excess area','FontSize',20,'FontWeight','bold','Color','k')
    pause(0.1)

    figure (8)
    title('Density vs position')
    plot(C(:,1),phi,'--bo','Markersize',4,'LineWidth',2)
    xlabel('R','FontSize',20,'FontWeight','bold','Color','k')
    ylabel('Density','FontSize',20,'FontWeight','bold','Color','k')

end

figure(990)
plot(time(2:end),max_phi_hist)
xlabel('time')
ylabel('maxphi')
figure(991)
plot(time(2:end),max_psi_hist)
xlabel('time')
ylabel('maxpsi')
figure(992)
plot(time(2:end),max_V_c_hist)
xlabel('time')
ylabel('maxVc')


%% osave =1 would save the output file while otherwise not
osave = 0;
output_file = 'Output_Files/Bud_Evolved_Branched_Actin.mat';

if osave == 1
    save(output_file);
end

return

