GPT tut interaction 20110530 first edit

Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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

Section 2
        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 10
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

Section 11
        % In this section you may modify the mesh in any way whatsoever.

Section 12
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 13
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 14
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 15
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );
Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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

Section 2
        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 10
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

Section 11
        % In this section you may modify the mesh in any way whatsoever.

Section 12
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 13
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 14
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 15
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );
Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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;

Section 2
        %%% USER CODE: INITIALISATION

        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


        %%% END OF USER CODE: INITIALISATION

Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
Section 10
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
Section 11
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 12
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

        % In this section you may modify the mesh in any way whatsoever.

Section 13
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 14
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 15
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 16
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );
Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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;

Section 2
        %%% USER CODE: INITIALISATION

        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


        %%% END OF USER CODE: INITIALISATION

Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
Section 10
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
Section 11
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 12
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

        % In this section you may modify the mesh in any way whatsoever.

Section 13
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 14
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 15
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 16
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );
Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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;

Section 2
        %%% USER CODE: INITIALISATION

        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


        %%% END OF USER CODE: INITIALISATION

Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
Section 10
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
Section 11
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 12
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

        % In this section you may modify the mesh in any way whatsoever.

Section 13
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 14
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 15
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 16
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );
Section 1
    function m = gpt_tut_interaction_20110530( m )
        %m = gpt_tut_interaction_20110530( m )
        %   Morphogen interaction function.
        %   Written at 2011-05-30 09:42:07.
        %   GFtbox revision 0, .

        % 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;

Section 2
        %%% USER CODE: INITIALISATION

        % In this section you may modify the mesh in any way whatsoever.
        if (Steps(m)==0) && m.globalDynamicProps.doinit % First iteration
            % Zero out a lot of stuff to create a blank slate.
            % If no morphogens are set in the GUI it may be useful to
            % zero some arrays by uncommenting the following.
            % m.morphogens(:) = 0;
            % m.morphogenclamp(:) = 0;
            % m.mgen_production(:) = 0;
            % m.mgen_absorption(:) = 0;
            % m.seams(:) = false;
            % m.mgen_dilution(:) = false;

            % Set up names for variant models.  Useful for running multiple models on a cluster.
            m.userdata.ranges.modelname.range = { 'MODEL1', 'MODEL2' };  % CLUSTER
            m.userdata.ranges.modelname.index = 1;                       % CLUSTER
        end
        modelname = m.userdata.ranges.modelname.range{m.userdata.ranges.modelname.index};  % CLUSTER
        disp(sprintf('\nRunning %s model %s\n',mfilename, modelname));
        switch modelname
            case 'MODEL1'
                % Set up the parameters (e.g. mutations) for this model here.
            case 'MODEL2'
                % Set up the parameters (e.g. mutations) for this model here.
            otherwise
                % If you reach here, you probably forgot a case.
        end

        % More examples of code for all iterations.

        % Set priorities for simultaneous plotting of multiple morphogens, if desired.
        % m = leaf_mgen_plotpriority( m, {'MGEN1', 'MGEN2'}, [1,2], [0.5,0.75] );

        % Set colour of polariser gradient arrows.
        % m = leaf_plotoptions(m,'highgradcolor',[0,0,0],'lowgradcolor',[1,0,0]);

        % setup a multiplot of the following morphogens
        % m = leaf_plotoptions( m, 'morphogen', {'V_PROFILE1','V_PROFILE2','KAPAR','S_LEFTRIGHT'});

        % to plot polariser on the A side and resultant areal growth rate on the B side:
        % m = leaf_plotoptions( m, 'morphogenA', 'POLARISER', ...
        %                      'outputquantityB', 'resultantgrowthrate', ...
        %                      'outputaxesB', 'areal' );

        % monitor properties of vertices must be done here - so that it reports newly equilibrated levels
        % m=leaf_profile_monitor(m,... % essential
        %         'REGIONLABELS',{'V_PROFILE1','V_PROFILE2'},... % essential
        %         'MORPHOGENS',{'S_LEFTRIGHT','S_CENTRE'},... % optional  (one element per REGIONLABEL)
        %         'VERTLABELS',false,'FigNum',1,'EXCEL',true,'MODELNAME',modelname); % optional (file in snapshots directory')


        %%% END OF USER CODE: INITIALISATION

Section 3
        %%% 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_plusorg_i,id_plusorg_p,id_plusorg_a,id_plusorg_l] = getMgenLevels( m, 'ID_PLUSORG' );
        [id_negorg_i,id_negorg_p,id_negorg_a,id_negorg_l] = getMgenLevels( m, 'ID_NEGORG' );
        [id_a_i,id_a_p,id_a_a,id_a_l] = getMgenLevels( m, 'ID_A' );
        [id_b_i,id_b_p,id_b_a,id_b_l] = getMgenLevels( m, 'ID_B' );

        % Mesh type: circle
        %          centre: 0
        %       circumpts: 48
        %       coneangle: 0
        %         dealign: 0
        %          height: 0
        %        innerpts: 0
        %      randomness: 0.1
        %           rings: 6
        %         version: 1
        %          xwidth: 0.2
        %          ywidth: 0.2

        %            Morphogen   Diffusion   Decay   Dilution   Mutant
        %            -------------------------------------------------
        %                KAPAR        ----    ----       ----     ----
        %                KAPER        ----    ----       ----     ----
        %                KBPAR        ----    ----       ----     ----
        %                KBPER        ----    ----       ----     ----
        %                 KNOR        ----    ----       ----     ----
        %            POLARISER        ----    ----       ----     ----
        %            STRAINRET        ----    ----       ----     ----
        %               ARREST        ----    ----       ----     ----
        %           ID_PLUSORG        ----    ----       ----     ----
        %            ID_NEGORG        ----    ----       ----     ----
        %                 ID_A        ----    ----       ----     ----
        %                 ID_B        ----    ----       ----     ----


        %%% USER CODE: MORPHOGEN INTERACTIONS

        % In this section you may modify the mesh in any way that does not
Section 4
        % alter the set of nodes.

        if (Steps(m)==0) && m.globalDynamicProps.doinit  % Initialisation code.
            % Put any code here that should only be performed at the start of
            % the simulation, for example, to set up initial morphogen values.
            % m.nodes is the set of vertex positions, an N by 3 array if there
            % are N vertices.  Row number K contains the X, Y, and Z
            % coordinates of the Kth vertex. To obtain a list of the X
            % coordinates of every vertex, write m.nodes(:,1).  The Y
            % coordinates are given by m.nodes(:,2) and the Z coordinates by
            % m.nodes(:,3).

            % Set up a morphogen promoter (_p suffix) region where x values are minimum
            % id_prox_p(m.nodes(:,1)==min(m.nodes(:,1)))=1;
            % if the morphogen level (_l suffix) is to be used in this iteration
            % set the level using the morphogen activity (_a suffix).
            % id_prox_l=id_prox_p * id_prox_a; % when a mutation is specified in the GUI
            % the activity (_a) is set to zero

            id_a_p(m.nodes(:,1)<-0.03)=1; % setup region for A where identity factor A is represented by id_a_p
            id_b_p(m.nodes(:,2)<-0.01)=1; % setup region for B

            % One way to set up a morphogen gradient is by ...
            % Setting up a gradient by clamping the ends (execute only once)
            % P=id_prox_p;
            % m.morphogenclamp( ((id_prox_p==1)|(id_dist_p==1)), polariser_i ) = 1;
            % m = leaf_mgen_conductivity( m, 'POLARISER', 0.01 );  %specifies the diffusion rate of polariser
            % m = leaf_mgen_absorption( m, 'POLARISER', 0.1 );     % specifies degradation rate of polariser

            % Fixing vertices, i.e. fix z for the base to prevent base from moving up or down
            % m=leaf_fix_vertex(m,'vertex',find(id_prox_p==1),'dfs','z');

            % To cut the mesh, set a temporary morphogen to 1 in places to cut
            % seams=zeros(size(P));
            % seams(indexes to places to cut)=1;
            % m=leaf_set_seams(m,seams);

        end

Section 5
        % Second way to generate a gradient
        % generating (+) and sinking (-) a diffusing signal (in this case polariser)
        % m.mgen_production( :, polariser_i ) = + 5*s_spur_p - P .* id_dist_p;

        % Monitor growth by scattering discs that deform over time (c.f. inducing biological clones)
        % (CARE - if the canvas is flat ensure that Plot:Hide Thickness is true,
        % because a quirk of the Matlab z-buffer means that they can get hidden by mistake)
        %    if (340>realtime-dt) && (340<realtime+dt) % discs to be added at realtime==340
        %        m = leaf_makesecondlayer( m, ...  % This function adds discs that represent transformed cells.
        %            'mode', 'each', ...  % Make discs randomly scattered over the canvas.
        %            'relarea', 1/16000, ...   % Each discs has area was 1/16000 of the initial area of the canvas.
        %            'probpervx', 'V_FLOWER', ... % induce discs over whole canvas (V_FLOWER is 1 over whole canvas)
        %            'numcells',4500,...%number of discs (that will become ellipses)
        %            'sides', 6, ...  % Each discs is approximated as a 6-sided regular polygon.
        %            'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %            'colorvariation',1,... % Each disc is a random colour
        %            'add', true );  % These discs are added to any discs existing already
        %    end

Section 6
        % Directives for creating latex representation directly from Matlab code
        % not fully implemented yet but will use @@ directives
        % @@at t
        % @@before t
        % @@after t
        % @@between t1 t2

Section 7
        %     % If you want to define different phases according to the absolute
        %     % time, create a morphogen for each phase and modulate
        %     % expressions using the morphogen
        %     % like.  For example:
        %     if (realtime < 10)  % first growth phase
        %         f_firstgrowth_p = 1;
        %     else
        %         f_firstgrowth_p = 0;
        %     end
        %     if (realtime >= 10) % second growth phase
        %         f_secondgrowth_p = 1;
        %     else
        %         f_secondgrowth_p = 0;
        %     end
        %
        %     % If you want one morphogen to affect others only during a certain
        %     % phase, write something like:
        %
        %     mgen_a_p = f_firstgrowth_p .* (various terms); % will zero except in firstgrowth

Section 8
        % Code common to all models.
        % @@PRN Polariser Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@GRN Gene Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.
        % @@KRN Growth Regulatory Network
        % Every equation to be formatted should end with an at-at Eqn N comment.

Section 9
        % Code for specific models.
        switch modelname
Section 10
            case 'MODEL1'  % @@model MODEL1
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
Section 11
            case 'MODEL2'  % @@model MODEL2
                % @@PRN Polariser Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % P(:) = ...  % @@ Eqn xx
                % @@GRN Gene Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % @@KRN Growth Regulatory Network
                % Every equation to be formatted should end with an at-at Eqn N comment.
                % kapar_p(:) = 0;  % @@ Eqn xx
                % kaper_p(:) = 0;  % @@ Eqn xx
                % kbpar_p(:) = 0;  % @@ Eqn xx
                % kbper_p(:) = 0;  % @@ Eqn xx
                % knor_p(:)  = 0;  % @@ Eqn xx
                kapar_p(:) = id_a_l .* inh(1,id_b_l); % growth rate
                kaper_p(:) = kapar_p; % isotropic growth
                kbpar_p(:) = kapar_p; % same on both sides of the sheet
                kbper_p(:) = kapar_p; % same
                knor_p(:)  = 0;       % thickness not growing
            otherwise
                % If this happens, maybe you forgot a model.
        end

Section 12
        %%% 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_plusorg_i) = id_plusorg_p;
        m.morphogens(:,id_negorg_i) = id_negorg_p;
        m.morphogens(:,id_a_i) = id_a_p;
        m.morphogens(:,id_b_i) = id_b_p;

        %%% USER CODE: FINALISATION

        % In this section you may modify the mesh in any way whatsoever.

Section 13
        % If needed force FE to subdivide (increase number FE's) here
        % if realtime==280+dt
        % m = leaf_subdivide( m, 'morphogen','id_vent',...
        %       'min',0.5,'max',1,...
        %       'mode','mid','levels','all');
        % end
        % Cut the mesh along the seams (see above)
        % if m.userdata.CutOpen==1
        %    m=leaf_dissect(m);
        %    m.userdata.CutOpen=2;
        %    Relax accumulated stresses slowly i.e. 0.95 to 0.999
        %    m = leaf_setproperty( m, 'freezing', 0.999 );
        % end

        %%% END OF USER CODE: FINALISATION

    end


    %%% USER CODE: SUBFUNCTIONS

Section 14
    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.020000 );
        %    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', 1.000000 );
        %    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', 0.000841 );
        %    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', 0.010000 );
        %    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, 'maxBioBcells', 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, 'freezing', 0.000000 );
        %    m = leaf_setproperty( m, 'canceldrift', false );
        %    m = leaf_setproperty( m, 'mgen_interaction',  );
        %    m = leaf_setproperty( m, 'mgen_interactionName', 'gpt_tut_interaction_20110530' );
        %    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', 'norm' );
        %    m = leaf_setproperty( m, 'diffusiontolerance', 0.000010 );
        %    m = leaf_setproperty( m, 'allowsparse', true );
        %    m = leaf_setproperty( m, 'maxIters', 0 );
        %    m = leaf_setproperty( m, 'maxsolvetime', 1000.000000 );
        %    m = leaf_setproperty( m, 'cgiters', 0 );
        %    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', false );
        %    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, 'jiggleProportion', 1.000000 );
        %    m = leaf_setproperty( m, 'cvtperiter', 0.200000 );
        %    m = leaf_setproperty( m, 'boingNeeded', false );
        %    m = leaf_setproperty( m, 'initialArea', 0.031326 );
        %    m = leaf_setproperty( m, 'bendunitlength', 0.176992 );
        %    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' );
        %    m = leaf_setproperty( m, 'modelname', 'GPT_tut_interaction_20110530' );
        %    m = leaf_setproperty( m, 'allowsave', true );
        %    m = leaf_setproperty( m, 'addedToPath', false );
        %    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', 0 );
        %    m = leaf_setproperty( m, 'coderevisiondate',  );
        %    m = leaf_setproperty( m, 'modelrevision', 0 );
        %    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', 0.200000 );
    end

Section 15
    % Here you may write any functions of your own, that you want to call from
    % the interaction function, but never need to call from outside it.
    % Remember that they do not have access to any variables except those
    % that you pass as parameters, and cannot change anything except by
    % returning new values as results.
    % Whichever section they are called from, they must respect the same
    % restrictions on what modifications they are allowed to make to the mesh.

    % For example:

Section 16
    % function m = do_something( m )
    %   % Change m in some way.
    % end

    % Call it from the main body of the interaction function like this:
    %       m = do_something( m );