Measuring using interaction function subroutines

Return to tutorials

Measuring growth of a model

Download a project illustrating measurements GPT_MeasureGrowth_20121203

The GPT_MeasureGrowth_20121203 interaction function is shown at the bottom of this page.

There are several ways to measure growth on the model. For example, there are two auxiliary functions:

[out1, out2, out3, out4] = leaf_computeGrowthFromDisplacements( m, displacements, time, ... )
   Compute the growth rate in each FE of m given the displacements of all
   the prism nodes and the time over which these displacements happened.
   This assumes that m is in the state after these
   displacements have happened.  The results returned depend on the
   options and on how many output arguments are given.  Note that unlike
   most leaf_* commands, m itself is not one of the output arguments,
   since this procedure does not need to modify m.

function  m=leaf_profile_monitor(m,realtime,RegionLabels,Morphogens,start_figno)
monitor morphogen levels at a set of vertices
m, mesh
RegionLabels, vertices to be monitored as designated by cell array of strings, i.e. region labels
Morphogens, cell array of strings, i.e. uppercase morphogen names to
   be monitored. There should be one RegionLabels string for each
   Morphogens string
Vertlabels, if true then display vertex numbers in each region label on
         the mesh default false
start_figno, default figure 1 (Must open a fig even if just monitoring to file)

MonData, optional output of data structure

• Measuring growth
• 

  (A) A single row of vertices (id_mid) labels a line along which measurements will be made. Here, it provides a source for s_mid.

• 

  (B) Kapar is promoted by s_a. Here, after 10 steps. We are going to view growth along id_mid.

• 

  (C) It has been through the interaction function 9 times. This shows the growth by the time it has reached the end of the 9th pass. Colour shows growth rate over the last time step and lines indicate the major axis.

• 

  (D) Spreadsheet showing distances along the line each time the system went through the interaction function. Line 3 shows the time. Column of numbers (lines 4 to 19) shows the distances of each vertex along i_mid. The spreadsheet is saved to the project snapshots directory.

• 

  (E) Left panel: distances between successive vertices along id_mid. Right panel: growth at each vertex along id_mid.

function m = gpt_measuregrowth_20121203( m )
%m = gpt_measuregrowth_20121203( m )
%   Morphogen interaction function.
%   Written at 2012-12-04 15:01:27.
%   GFtbox revision 4418, .

% The user may edit any part of this function between delimiters
% of the form "USER CODE..." and "END OF USER CODE...".  The
% delimiters themselves must not be moved, edited, deleted, or added.

    if isempty(m), return; end

    fprintf( 1, '%s found in %s\n', mfilename(), which(mfilename()) );

    try
        m = local_setproperties( m );
    catch
    end

    realtime = m.globalDynamicProps.currenttime;


%%% USER CODE: INITIALISATION

%%% END OF USER CODE: INITIALISATION


%%% SECTION 1: ACCESSING MORPHOGENS AND TIME.
%%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT.

    if isempty(m), return; end

    setGlobals();
    global gNEW_KA_PAR gNEW_KA_PER gNEW_KB_PAR gNEW_KB_PER
    global gNEW_K_NOR gNEW_POLARISER gNEW_STRAINRET gNEW_ARREST
    dt = m.globalProps.timestep;
    polariser_i = gNEW_POLARISER;
    P = m.morphogens(:,polariser_i);
    [kapar_i,kapar_p,kapar_a,kapar_l] = getMgenLevels( m, 'KAPAR' );
    [kaper_i,kaper_p,kaper_a,kaper_l] = getMgenLevels( m, 'KAPER' );
    [kbpar_i,kbpar_p,kbpar_a,kbpar_l] = getMgenLevels( m, 'KBPAR' );
    [kbper_i,kbper_p,kbper_a,kbper_l] = getMgenLevels( m, 'KBPER' );
    [knor_i,knor_p,knor_a,knor_l] = getMgenLevels( m, 'KNOR' );
    [strainret_i,strainret_p,strainret_a,strainret_l] = getMgenLevels( m, 'STRAINRET' );
    [arrest_i,arrest_p,arrest_a,arrest_l] = getMgenLevels( m, 'ARREST' );
    [id_mid_i,id_mid_p,id_mid_a,id_mid_l] = getMgenLevels( m, 'ID_MID' );
    [s_mid_i,s_mid_p,s_mid_a,s_mid_l] = getMgenLevels( m, 'S_MID' );

% Mesh type: rectangle
%            base: 0
%          centre: 0
%      randomness: 0.1
%         version: 1
%           xdivs: 16
%          xwidth: 2
%           ydivs: 16
%          ywidth: 2

%            Morphogen    Diffusion   Decay   Dilution   Mutant
%            --------------------------------------------------
%                KAPAR         ----    ----       ----     ----
%                KAPER         ----    ----       ----     ----
%                KBPAR         ----    ----       ----     ----
%                KBPER         ----    ----       ----     ----
%                 KNOR         ----    ----       ----     ----
%            POLARISER         ----    ----       ----     ----
%            STRAINRET         ----    ----       ----     ----
%               ARREST         ----    ----       ----     ----
%               ID_MID         ----    ----       ----     ----
%                S_MID         0.01    ----       ----     ----


%%% USER CODE: MORPHOGEN INTERACTIONS

    if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
        % assign a morphogen down the middle of the y axis
        id_mid_p((m.nodes(:,2)<0.1)&(m.nodes(:,2)>-0.1))=1;
        s_mid_p(:)=0;
        s_mid_p=id_mid_p;
        m.morphogenclamp( (id_mid_p==1), s_mid_i ) = 1;
        m = leaf_mgen_conductivity( m, 's_mid', 0.01 );  %specifies the diffusion rate of polariser    
        m = leaf_mgen_absorption( m, 's_mid', 0.0 );     % specifies degradation rate of polariser
    end
        kapar_p=15 * pro(1,s_mid_p);
        kbpar_p=15 * pro(1,s_mid_p);
        kaper_p=10 * pro(1,s_mid_p);
        kbper_p=10 * pro(1,s_mid_p);
        knor_p =0;
    % calculate growth over dt interval
    if isfield(m.userdata,'oldpos')
        displacements = m.prismnodes - m.userdata.oldpos;
        [growth,gf] = leaf_computeGrowthFromDisplacements( m, displacements, ...
            realtime - m.userdata.starttime,'axisorder', 'maxminnor', ...
            'anisotropythreshold', 0.05);
        % plot resultant areal growth rates over dt (i.e. one step).
        m = leaf_plotoptions( m, 'pervertex',perFEtoperVertex(m,sum(growth(:,1:2),2)) ,'perelementaxes', gf(:,1,:), 'drawtensoraxes', true );
        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
    end
    m.userdata.oldpos = m.prismnodes;
    m.userdata.starttime = realtime;
    % monitor growth in a Figure and write excel file
    [m,MonData]=leaf_profile_monitor(m,... % essential
        'REGIONLABELS',{'ID_MID'},... % essential
        'VERTLABELS',false,'FigNum',1,'EXCEL',true); % optional (file in snapshots directory')
    
% if sum(ismember(output,realtime)) && m.userdata.output ==1
%     % Output - plot an image at high resolution   
%     path = fileparts(which(m.globalProps.modelname));
%     [m,ok] = leaf_snapshot( m,[path,filesep,'snapshots',filesep,modelname,'_',modeltype,'.png'], 'resolution',[]); %
% 
% end

%%% END OF USER CODE: MORPHOGEN INTERACTIONS

%%% SECTION 3: INSTALLING MODIFIED VALUES BACK INTO MESH STRUCTURE
%%% AUTOMATICALLY GENERATED CODE: DO NOT EDIT.
    m.morphogens(:,polariser_i) = P;
    m.morphogens(:,kapar_i) = kapar_p;
    m.morphogens(:,kaper_i) = kaper_p;
    m.morphogens(:,kbpar_i) = kbpar_p;
    m.morphogens(:,kbper_i) = kbper_p;
    m.morphogens(:,knor_i) = knor_p;
    m.morphogens(:,strainret_i) = strainret_p;
    m.morphogens(:,arrest_i) = arrest_p;
    m.morphogens(:,id_mid_i) = id_mid_p;
    m.morphogens(:,s_mid_i) = s_mid_p;

%%% USER CODE: FINALISATION

%%% END OF USER CODE: FINALISATION

end


%%% USER CODE: SUBFUNCTIONS
%%% USER CODE: SUBFUNCTIONS


function m = local_setproperties( m )
% This function is called at time zero in the INITIALISATION section of the
% interaction function.  It provides commands to set each of the properties
% that are contained in m.globalProps.  Uncomment whichever ones you would
% like to set yourself, and put in whatever value you want.
%
% Some of these properties are for internal use only and should never be
% set by the user.  At some point these will be moved into a different
% component of m, but for the present, just don't change anything unless
% you know what it is you're changing.

%    m = leaf_setproperty( m, 'trinodesvalid', true );
%    m = leaf_setproperty( m, 'prismnodesvalid', true );
%    m = leaf_setproperty( m, 'thicknessRelative', 0.200000 );
%    m = leaf_setproperty( m, 'thicknessArea', 0.000000 );
%    m = leaf_setproperty( m, 'thicknessMode', 'physical' );
%    m = leaf_setproperty( m, 'activeGrowth', 1.000000 );
%    m = leaf_setproperty( m, 'displayedGrowth', 10 );
%    m = leaf_setproperty( m, 'displayedMulti', [] );
%    m = leaf_setproperty( m, 'allowNegativeGrowth', true );
%    m = leaf_setproperty( m, 'usePrevDispAsEstimate', true );
%    m = leaf_setproperty( m, 'perturbInitGrowthEstimate', 0.000010 );
%    m = leaf_setproperty( m, 'perturbRelGrowthEstimate', 0.010000 );
%    m = leaf_setproperty( m, 'perturbDiffusionEstimate', 0.000100 );
%    m = leaf_setproperty( m, 'resetRand', false );
%    m = leaf_setproperty( m, 'mingradient', 0.000000 );
%    m = leaf_setproperty( m, 'relativepolgrad', false );
%    m = leaf_setproperty( m, 'usefrozengradient', true );
%    m = leaf_setproperty( m, 'userpolarisation', false );
%    m = leaf_setproperty( m, 'thresholdsq', 1.000000 );
%    m = leaf_setproperty( m, 'splitmargin', 1.400000 );
%    m = leaf_setproperty( m, 'splitmorphogen',  );
%    m = leaf_setproperty( m, 'thresholdmgen', 0.500000 );
%    m = leaf_setproperty( m, 'bulkmodulus', 1.000000 );
%    m = leaf_setproperty( m, 'unitbulkmodulus', true );
%    m = leaf_setproperty( m, 'poissonsRatio', 0.300000 );
%    m = leaf_setproperty( m, 'starttime', 0.000000 );
%    m = leaf_setproperty( m, 'timestep', 1.000000 );
%    m = leaf_setproperty( m, 'timeunitname',  );
%    m = leaf_setproperty( m, 'distunitname', 'mm' );
%    m = leaf_setproperty( m, 'scalebarvalue', 0.000000 );
%    m = leaf_setproperty( m, 'validateMesh', true );
%    m = leaf_setproperty( m, 'rectifyverticals', false );
%    m = leaf_setproperty( m, 'allowSplitLongFEM', true );
%    m = leaf_setproperty( m, 'longSplitThresholdPower', 0.000000 );
%    m = leaf_setproperty( m, 'allowSplitBentFEM', false );
%    m = leaf_setproperty( m, 'allowSplitBio', true );
%    m = leaf_setproperty( m, 'allowFlipEdges', false );
%    m = leaf_setproperty( m, 'allowElideEdges', true );
%    m = leaf_setproperty( m, 'mincellangle', 0.200000 );
%    m = leaf_setproperty( m, 'alwaysFlat', 0.000000 );
%    m = leaf_setproperty( m, 'flattenforceconvex', true );
%    m = leaf_setproperty( m, 'flatten', false );
%    m = leaf_setproperty( m, 'flattenratio', 1.000000 );
%    m = leaf_setproperty( m, 'useGrowthTensors', false );
%    m = leaf_setproperty( m, 'plasticGrowth', false );
%    m = leaf_setproperty( m, 'totalinternalrotation', 0.000000 );
%    m = leaf_setproperty( m, 'stepinternalrotation', 2.000000 );
%    m = leaf_setproperty( m, 'showinternalrotation', false );
%    m = leaf_setproperty( m, 'performinternalrotation', false );
%    m = leaf_setproperty( m, 'internallyrotated', false );
%    m = leaf_setproperty( m, 'maxFEcells', 0 );
%    m = leaf_setproperty( m, 'inittotalcells', 0 );
%    m = leaf_setproperty( m, 'bioApresplitproc',  );
%    m = leaf_setproperty( m, 'bioApostsplitproc',  );
%    m = leaf_setproperty( m, 'maxBioAcells', 0 );
%    m = leaf_setproperty( m, 'colors', (6 values) );
%    m = leaf_setproperty( m, 'colorvariation', 0.050000 );
%    m = leaf_setproperty( m, 'colorparams', (12 values) );
%    m = leaf_setproperty( m, 'biocolormode', 'auto' );
%    m = leaf_setproperty( m, 'freezing', 0.999000 );
%    m = leaf_setproperty( m, 'canceldrift', false );
%    m = leaf_setproperty( m, 'mgen_interaction', (unknown type function_handle) );
%    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_demosubdivision_20121114' );
%    m = leaf_setproperty( m, 'allowInteraction', true );
%    m = leaf_setproperty( m, 'interactionValid', true );
%    m = leaf_setproperty( m, 'gaussInfo', (unknown type struct) );
%    m = leaf_setproperty( m, 'stitchDFs', [] );
%    m = leaf_setproperty( m, 'D', (36 values) );
%    m = leaf_setproperty( m, 'C', (36 values) );
%    m = leaf_setproperty( m, 'G', (6 values) );
%    m = leaf_setproperty( m, 'solver', 'cgs' );
%    m = leaf_setproperty( m, 'solverprecision', 'double' );
%    m = leaf_setproperty( m, 'solvertolerance', 0.001000 );
%    m = leaf_setproperty( m, 'solvertolerancemethod', 'max' );
%    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
%    m = leaf_setproperty( m, 'allowsparse', true );
%    m = leaf_setproperty( m, 'maxIters', 3.000000 );
%    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
%    m = leaf_setproperty( m, 'cgiters', 0.000000 );
%    m = leaf_setproperty( m, 'simsteps', 0 );
%    m = leaf_setproperty( m, 'stepsperrender', 0 );
%    m = leaf_setproperty( m, 'growthEnabled', true );
%    m = leaf_setproperty( m, 'diffusionEnabled', true );
%    m = leaf_setproperty( m, 'flashmovie', false );
%    m = leaf_setproperty( m, 'makemovie', 0.000000 );
%    m = leaf_setproperty( m, 'moviefile',  );
%    m = leaf_setproperty( m, 'codec', 'None' );
%    m = leaf_setproperty( m, 'autonamemovie', true );
%    m = leaf_setproperty( m, 'overwritemovie', false );
%    m = leaf_setproperty( m, 'framesize', [] );
%    m = leaf_setproperty( m, 'mov', [] );
%    m = leaf_setproperty( m, 'boingNeeded', false );
%    m = leaf_setproperty( m, 'initialArea', 4.000000 );
%    m = leaf_setproperty( m, 'bendunitlength', 2.000000 );
%    m = leaf_setproperty( m, 'targetRelArea', 1.000000 );
%    m = leaf_setproperty( m, 'defaultinterp', 'min' );
%    m = leaf_setproperty( m, 'readonly', false );
%    m = leaf_setproperty( m, 'projectdir', 'D:\ab\Matlab stuff\Growth models' );
%    m = leaf_setproperty( m, 'modelname', 'GPT_DemoSubdivision_20121114' );
%    m = leaf_setproperty( m, 'allowsave', 1.000000 );
%    m = leaf_setproperty( m, 'addedToPath', true );
%    m = leaf_setproperty( m, 'bendsplit', 0.300000 );
%    m = leaf_setproperty( m, 'usepolfreezebc', false );
%    m = leaf_setproperty( m, 'dorsaltop', true );
%    m = leaf_setproperty( m, 'defaultazimuth', -45.000000 );
%    m = leaf_setproperty( m, 'defaultelevation', 33.750000 );
%    m = leaf_setproperty( m, 'defaultroll', 0.000000 );
%    m = leaf_setproperty( m, 'defaultViewParams', (unknown type struct) );
%    m = leaf_setproperty( m, 'comment',  );
%    m = leaf_setproperty( m, 'legendTemplate', '%T: %q\n%m' );
%    m = leaf_setproperty( m, 'bioAsplitcells', true );
%    m = leaf_setproperty( m, 'bioApullin', 0.142857 );
%    m = leaf_setproperty( m, 'bioAfakepull', 0.202073 );
%    m = leaf_setproperty( m, 'interactive', false );
%    m = leaf_setproperty( m, 'coderevision', 4351.000000 );
%    m = leaf_setproperty( m, 'coderevisiondate',  );
%    m = leaf_setproperty( m, 'modelrevision', 0.000000 );
%    m = leaf_setproperty( m, 'modelrevisiondate',  );
%    m = leaf_setproperty( m, 'savedrunname',  );
%    m = leaf_setproperty( m, 'savedrundesc',  );
%    m = leaf_setproperty( m, 'vxgrad', (108 values) );
%    m = leaf_setproperty( m, 'lengthscale', 2.000000 );
end