% 2D VERTEX MODEL OF DROSOPHILA VENTRAL FURROW FORMATION
% MATLAB version with utilization of mex compilers
% By Oleg Polyakov | Princeton University | 2008 - 2013

% This model was developed to study the process of ventral furrow formation
% in Drosophila Melanogaster.

% In this model we consider the consider the midplane cross section of the
% emrbyo at the onset of gasrtulation. This corresponds to a geometry of 80
% cells connected in a ring. The volumes of the yolk, the cellular
% volumes, as well as the apical, basal, and lateral sides of all the cells
% have an "elastic" energy equation of the form K*(X - Xo)^2. 

% The input energy comes entirely from apical constriction. It is a
% Gaussian distribution centered on the ventral surface with a width of
% sigma.

% The supporting mex files are the following:
% INIT_C : initializes the cross sectional geometry, the equilibrium values
% of the cells and apical constriction

% CALC_GRAD_C: calculates the gradient by taking the derivative of the
% energy equation

% CALC_GRAD_C_N_STEPS: calculates the gradient and moves the verticies
% along the gradient with a magnitude of P.beta. This is repeated P.N_steps
% times. If equilibrium is reached then the flag P.Done is set to 0.

% GEOM_OUT_C: outputs the current geometry parameters for each cell.


%% Init
clearvars -except ID_KL ID_KB;

% Flag to see when simulation comes to a stop
P.Done = 1;

% File to save the data
savefile_path = [pwd, '/'];
savefile = [savefile_path, 'test.mat'];

% Outside and inner radius of the emrbyo
P.R_out = 14; P.R_in = 8.9;

% Apical constriction
P.Alpha = [1, 5000, 0]; % [linear, quandratic, basal contriction (not used)]
P.S = [4, 3, 0]; % sigma for above (zero means flat)

%KL_val = KL_vector(ID_KL);
KL_val = 20.0;
KA_val = 5.0;
KB_val = 1.0;
KV_val = 5000;
KY_val = [1, 0];

% Simulation runs
P.N_steps = 100000; P.beta = .0025/KV_val;P.N_steps_end = 10000-1;

% Inititialize values
[Ax, Ay, Bx, By, P.Ao_con, P.Ao, P.KA_lin, P.KA_con, P.KA, P.Bo, P.KB_lin, P.KB,...
    P.Lo, P.KL, P.Vo, P.KV, P.Yo, P.KY] = INIT_C(P.R_out, P.R_in, P.Alpha, P.S,...
    KA_val, KB_val, KL_val, KV_val, KY_val);

n_meso_Ao_drop = 12;K_meso_Ao_drop = 7;
P.Ao_con = P.Ao_con.*(1 - 1./(1 + (((1:80)-40.5)./K_meso_Ao_drop).^n_meso_Ao_drop));

%n_ecto_Lo = 6; K_ecto_Lo = 10;Lo_factor = .75;
%P.Lo = P.Lo.*((1-Lo_factor)*1./(1 + (((1:80)-41)./K_ecto_Lo).^n_ecto_Lo) + Lo_factor);

%% Init Gradients
P.count = 0; P.loop_count = 0;
[Gax, Gay, Gbx, Gby, E_old] = CALC_GRAD_C(Ax, Ay, Bx, By,...
    P.KA_lin, P.KA_con, P.Ao_con, P.KA, P.Ao, P.KL, P.Lo, P.KB, P.KB_lin, P.Bo, P.KV, P.Vo,...
    P.KY, P.Yo, P.R_out, P.R_out, P.beta); 

P.count = P.count + 1; 
[P.A(P.count, :), P.B(P.count, :), P.L(P.count, :), P.V(P.count, :), P.Y(P.count, :)] = GEOM_OUT_C(Ax, Ay, Bx, By);
Vxm{P.count} = [Ax;Bx]; Vym{P.count} = [Ay;By]; Gxm{P.count} = [Gax;Gbx]; Gym{P.count} = [Gay;Gby];    

% Run simulation until equilibrium is reached
while (P.Done),
    [Gax, Gay, Gbx, Gby, P.Done] = CALC_GRAD_C_N_STEPS(Ax, Ay, Bx, By,...
        P.KA_lin, P.KA_con, P.Ao_con, P.KA, P.Ao, P.KL, P.Lo, P.KB, P.KB_lin, P.Bo, P.KV, P.Vo,...
        P.KY, P.Yo, P.R_out, P.R_out, P.beta, P.N_steps, E_old); P.loop_count = P.loop_count + 1;
    
    disp(['Elongation is ', num2str(P.L(end, 40)/5.1), '  Volume change is ', num2str(min(P.V(end, :)/P.Vo(1)))]);
    P.count = P.count + 1;
    [P.A(P.count, :), P.B(P.count, :), P.L(P.count, :), P.V(P.count, :), P.Y(P.count, :)] = GEOM_OUT_C(Ax, Ay, Bx, By);
    Vxm{P.count} = [Ax;Bx]; Vym{P.count} = [Ay;By]; Gxm{P.count} = [Gax;Gbx]; Gym{P.count} = [Gay;Gby];
    
    if P.loop_count == P.N_steps_end, break; end
end

P.t_stop = P.loop_count;
%%{
if exist(savefile, 'file'), delete(savefile); end
save(savefile, 'Vxm', 'Vym', 'Gxm', 'Gym', 'P');

%%}


