Rubio_etal_2021/Code/Scripts/MAIN_Axon_Model_Pepo_1.m

clear all

%% Main Script to Generate Axons

%  Copyright (C) INRIA
%  This file is released under the CeCILL-B license. See LICENSE.txt.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%This Script allows to automatically perfom simulations changing different parameters. 

%The main file 'FileName' will be saved for each simulation.

%FileName contains the cell-type variable Axons will all the information regarding the simulations.

%Mutants can be generated within the population that are: slower (changing mut_slow), that have different parmeter values (Numel(Alpha) and Numel(Beta)=2) and that present different probabilities of mechanical branching (IMP_ran(2)). Other possibilities may be included in the code.

%Axon cell structure: 

    %-each neurite (main axon or branch) is described by a structure which represents anelement of the cell-type variable Axons. 
    %-each neurite is structured as follows:
    
%            mut: 0                             %0 means WT, 1 means mutant
%            level: 0                           %0 means main axon, 1 means first order branch and 2 second order branch.
%               bp: 0                           %if the neurite is a branch, is the branching point. If not is 0.
%               ND: [17x3 double]               %All the branching points
%           ND_dyn: [13x3 double]               %Type II (short and dynamic) Branching points
%               DB: {1x19 cell}                 %cell-type variable. Each element has the information of all the Type II branches that where created. The boolean 'stable' indicates if the branch has been stabilized or not. The variable 'origin' indicates ifit was stabilized by the contact with another branch tip (origin='BT') or branching point (origin='BP')
%     branch_index: [19x1 double]               %Element of the Axon cell-type variable where each branch is saved
%             time: 63                          %Time when theneurite reached the end
%       root_index: -1                          % -1 if main axon, if branch is the element of the Axon cell-type variable where each branch is saved
%       Parameters: [7.4450 1.6650]             % Parameters(1)=Alpha and Parameters(2)=Beta.
%             axon: [132x3 double]              %3D neurite path
%               th: [132x2 double]              %Theta values: th(:,1) for the XY plane and th(:,2) for the Z direction
%        SpaceTime: [63x3 double]               %size(SpaceTime,1) depicts the time points while the neurite grew. For each time point, the 3D vector depicts the spatiallocation of the axon tip. This allows to reconstruct the path in time.
%          counter: 84                          %final counter value of the neurite


%Code generated by Agustina Razetti (Université Côte d’Azur, Inria, CNRS, I3S, France). Contact: arazetti@unice.fr
%For any publication using this code, cite: Modelling axon growth fom in vivo data reveals the importance of (..)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NN=1; %Number of simulations -> many simulations can be performed. Identical or varying one or more parameters each time.

tic
for i=1:NN
    
%     close all

    Num_ax = 1;%Number of axons.
    
    % -------------  Parámetros por defecto  ---------------

    mut_slow=[];%if empty: nothing changes. if not: mut_slow is the equivalent to nmax_nr for mutants (instead of speed [nmax nr]=[6 2], axons may be slower, for example [3 1]
    Diameter=0.4;% axon diameter. (0.23 for gamma)
    Param_Error=[];%Param_Error=[]; means fixed values of Alpha and Beta. If Param_Error=[a b,c d], the parameters of each neuron are thrown from a Bivariate Gaussian of mean(Alpha Beta) and covariance matrix Param_Error.
    %Param_Error(:,:,2)=[]; %If wild types and mutants have different values, Param_Error(:,:,1) is the covariance matrix for wild type neurons and Param_Error(:,:,2) for mutants
    IMP_ran=[1 1]; %Probability of branching upon interaction (from 0-1). IMP_ran(1) for WT (typically equal to 1) and IMP_ran(2) for mutants (for example, equal to 0 for imp)
            % if Pb(1)>0 and IMP_ran(1)>0 both types of branches will appear: random and upon interaction
    % Lambda_b (Poisson number for branch distance) is set in Line 8 of code file 'Curve_Master_Branches_no_video_17'
    % Initial branch angle distribution can be set in Fio in line 274 of code file 'Curve_Master_Branches_no_video_17'
    num_mutants=0;%Number of mutants
%     FileName='try_wt';
    nexp=i; %experiment_counter (automatic)
    max_counter=140;%Counter number. Max. number of times a neurite may encounter a mechanical constraint. 


% ***********   PYRAMIDAL NEURONS (paper)   *****************

%     % -------------  PARÁMETROS CENTRALES MAIN AXON  ---------------
%    
%     Alpha = 10;%if numel(Alpha)=1 : every axon has the same value. if numel(Alpha)=2 : wild types will have Alpha(1) and mutants Alpha(2).
%     Beta = 30;%if numel(Beta)=1 : every axon has the same value. if numel(Beta)=2 : wild types will have Beta(1) and mutants Beta(2).
%     Pb = [0.3 0];%Pb(1): Probability of random branching (from 0-1). 
%              %Pb(2): by default equal to zero. If 1, the axon will not only try to branch with probabilty Pb(1) but will also do not advance during the corrspondent timepoint. It simulates random pausing.
    
% -------------  PARÁMETROS CENTRALES BRANCHES  ---------------
   
    Alpha = 50;%if numel(Alpha)=1 : every axon has the same value. if numel(Alpha)=2 : wild types will have Alpha(1) and mutants Alpha(2).
    Beta = 0.1;%if numel(Beta)=1 : every axon has the same value. if numel(Beta)=2 : wild types will have Beta(1) and mutants Beta(2).
    Pb = [0.2 0];%Pb(1): Probability of random branching (from 0-1). 
             %Pb(2): by default equal to zero. If 1, the axon will not only try to branch with probabilty Pb(1) but will also do not advance during the corrspondent timepoint. It simulates random pausing.
    
    % ATENCIÓN! si queremos modelizar las pyramidal neurons del paper necesitamos que Pb dependa de la x. 
    % Esto se controla con el parámetro 'Cond_x_Prob_Branching', en el script "Curve_master_branches_no_video_17.m":
    
    % -------------  Restricciones físicas: espacio disponible, distancia máxima, paso y nivel de ramificación ---------------
    
    Xi = 1;% Growth cavity. If Xi=1: a Tube is considered (Diameter set in line 14 of code file 'Curve_Simulation_branches_no_video_17'. If Xi=0: The Gamma axons medial lobe. 
    Xmax = 80; %Max. Travelled distance
    Drho = 1/3; %StepSize
    b_l = 2; %Max. branch order/level
    nmax_nr=[10 1];%nmax_nr(1): max. number of steps per time point (nmax). nmax_nr(2)=number of steps the neurite retracts upon mechanical interaction (nr)
    lambda_b = 3; % lambda_b es el parámetro 'MeanLambda' en "Curve_master_branches_no_video_17.m", línea 12.
    
    % -------------  Nombre del archivo ---------------
    
    FileName='Prueba_Pepo_1';

    % -------------  Definición del entorno ---------------
    
%     GRADIENTfile='GradientMap_Pepo_Pyramidal';% External attractive field to simulate the pyramidal neurons in the paper. 
    GRADIENTfile='GradMap_zero_for_TUBE_long';%Gradient Field (constant and equal zero.
%     GRADIENTfile='GradientMap_03_06'; %To simulate Gamma axons use : 'GradientMap_03_06.mat';
    
% ******************************************************************

                 
    % -------------  SIMULACIÓN ---------------

    SIMULATIONS_17_PEPO(Num_ax,Diameter,nmax_nr,FileName,nexp,IMP_ran,num_mutants,max_counter,mut_slow,[Alpha Beta],GRADIENTfile,Xi,Param_Error,Pb,b_l,Drho,Xmax, lambda_b);
         
end
toc

% %%  Plotting

clc
load(['exp_ 1 _ ', FileName, '.mat']);

figure
hold on
for k = 1:length(Axons)
    axon = Axons{1, k};
    if length(axon.axon) >= 10
        plot3(axon.axon(:, 1), axon.axon(:, 2), axon.axon(:, 3), 'linewidth', 1.5);
    end
end
xlabel('x');
ylabel('y');
zlabel('z');
title({['alpha = ', num2str(Alpha), '.  beta = ', num2str(Beta), '.  Pb = ', num2str(Pb(1)), '.  lambda b = ', num2str(lambda_b)], ['Xmax = ', num2str(Xmax), '.  Drho = ', num2str(Drho), '.  bl = ', num2str(b_l), '.  nmax = ', num2str(nmax_nr(1))]}, 'fontsize', 15);
axis equal