

disp(sprintf('mutant state cyc=%d dich=%d div=%d cenorg=%d dichdistor=%d',...
    id_cyc_a,id_dich_a,id_div_a,id_cenorg_a,v_dichdistorg_a))
OLD_DIFFUSION = false;
HOURS_PER_PLASTOCHRON = 10;
SPACERESCALE = 4;
pd_offset=0;

if realtime<340
    id_late_p(:)=0;
    id_early_p(:)=1;
else
    id_late_p(:)=1;
    id_early_p(:)=0;
end

%INITIALISATION
if realtime==22*HOURS_PER_PLASTOCHRON %23*HOURS_PER_PLASTOCHRON
    m.userdata.all_regions_ready=false; % controls dynamic subdivision
    m.userdata.donesubdivisions=false;
    m.userdata.extra_regions_ready=false;
    fprintf( 1, 'Zeroing morphogens at step %d, time %.3f.\n', ...
        Steps(m), realtime );
    m.morphogens=zeros(size(m.morphogens));
    m.mgen_production(:)=0;
    fprintf( 1, 'Initialising morphogens at step %d, time %.3f.\n', ...
        Steps(m), realtime );
    
    % First assign factors to regions (vertices) in the finite elements
    % mesh. These are shown in Figures 5 and 6.
    %
    % v_theta is the angle around the z axis.  This is for visualisation
    % purposes and does not model a biological morphogen.
    v_theta_p = atan2( m.nodes(:,2), m.nodes(:,1) )*(0.5/pi);
    
    % Find all the nodes at the base of the initial flower shape,
    % i.e. those whose z
    % coordinate is less than -0.5/SPACERESCALE.  Set the id_prox morphogen to be 1
    % there.  Since id_prox does not diffuse and is not absorbed, the
    % values set here persist unchanged during development.
    basenodes = m.nodes(:,3) < -0.1801;
    
    % Find the inner base nodes: those whose radius in the XY plane is
    % less than the maximum.
    xy = m.nodes(:,[1 2]);
    radiisq = sum(xy.*xy,2);
    maxrsq = max(radiisq);
    innerbasenodes = radiisq < 0.20*0.99;
    m.userdata.innerbasenodes=innerbasenodes;
    
    id_prox_p(basenodes) = 1;
    id_prox_l = id_prox_p;
    m = leaf_fix_vertex( m, 'vertex', basenodes, 'dfs', 'z' );
    
    % create region to be marked with clones
    v_flower_p=ones(size(v_flower_p));
    v_flower_p(basenodes)=0;
    
    % Find all the border nodes.
    borderedges = find( m.edgecells(:,2)==0 );
    bordernodes = unique( m.edgeends( borderedges, : ) );
    bordercells = m.edgecells( borderedges, 1 );
    thickbordernodes = unique( m.tricellvxs( bordercells, : ) );
    
    % Set id_dist to 1 at every node on the distal edge.
    id_dist_p( bordernodes ) = 1;
    id_dist_l = id_dist_p .* id_dist_a;
    
    
    maxlatheight = 0.15;
    latrimnodes = intersect( find(m.nodes(:,3) < maxlatheight), bordernodes );
    id_lat_p(:) = 0;
    id_lat_p(latrimnodes) = 1;
    petalborderv_theta = (1:5)*0.2 - 0.5;
    petalmidv_theta = petalborderv_theta - 0.1;
    v_thetatol = 0.002;
    petalbordernodes = false( size(m.nodes,1), 1 );
    petalmidnodes = false( size(m.nodes,1), 1 );
    for i=1:length(v_theta_p)
        th = v_theta_p(i);
        bth = abs(petalborderv_theta-th);
        if any( (bth < v_thetatol) | (1-bth < v_thetatol) )
            petalbordernodes(i) = true;
        end
        mth = abs(petalmidv_theta-th);
        if any( (mth < v_thetatol) | (1-mth < v_thetatol) )
            petalmidnodes(i) = true;
        end
    end
    id_lat_p(petalbordernodes) = 1;
    id_lat_p(innerbasenodes) = 0;
    id_med_p(:) = zeros(size(id_med_p));
    id_med_p(petalmidnodes) = 1;
    id_med_p(innerbasenodes) = 0;
    id_lat_l = id_lat_p .* id_lat_a;
    id_med_l = id_med_p .* id_med_a;
    
    id_secvein_p = 0.5 + 0.5*sin(20*pi*v_theta_p - pi/2);
    id_secvein_p = id_secvein_p > 0.80; % Thresholding to keep the region narrow.  If using a 0.5 of threshold continuous in P-D axis but very thick.
    %As the function of secvein will be to position the region of high Kper, let's concentrate on getting a narrow strip of secvein in the middle of each petal half
    id_secvein_l = id_secvein_p .* id_secvein_a;
    
    lobenodes = m.nodes(:,3) > -0.05;
    id_lobe_p(lobenodes) = 1;
    % the levels of clamped nodes (vertices) cannot change during
    % simulation
    m.morphogenclamp(lobenodes,id_lobe_i) = 1;
    
    id_tubenodes = m.nodes(:,3) <-0.03 & m.nodes(:,3) > - 0.18;
    id_tube_p(id_tubenodes) = 1;
    m.morphogenclamp(id_tubenodes,id_tube_i) = 1;
    
    id_lipnodes = m.nodes(:,3) > pd_offset-0.11 & m.nodes(:,3) < pd_offset+0.02;
    id_lip_p(id_lipnodes) = 1;
    m.morphogenclamp(id_lipnodes,id_lip_i) = 1;
    
    id_rimnodes = m.nodes(:,3) > pd_offset-0.07 & m.nodes(:,3) < pd_offset-0.03;
    id_rim_p(id_rimnodes) = 1;
    m.morphogenclamp(id_rimnodes,id_rim_i) = 1;
    
    % activate id_dich in dorsal domain
    id_dich_p=zeros(size(id_dich_p));
    id_dichnodes = (m.nodes(:,1) < -1.70/SPACERESCALE) & ~m.userdata.innerbasenodes;
    id_dich_p(id_dichnodes) = 1;
    id_dich_l=id_dich_p.*id_dich_a;
    
    % activate id_cyc in dorsal domain
    id_cyc_p=zeros(size(id_cyc_p));
    id_cycnodes = (m.nodes(:,1) < -0.153) & ~m.userdata.innerbasenodes;
    id_cyc_p(id_cycnodes) = 1;
    id_cyc_l=id_cyc_p.*id_cyc_a;
    
    f_seppetals_p(id_tube_p>0)=1; % all tube
    f_seppetals_p(id_lip_p>0)=1; % all lip
    f_seppetals_p=f_seppetals_p.*id_lat_p; % separate petals
    id_rimnodes = m.nodes(:,3) > pd_offset-0.01 & m.nodes(:,3) < pd_offset+0.01; % cut off lobes
    f_seplobes1_p(id_rimnodes)=1; % separate lobes from tube
    id_rimnodes = m.nodes(:,3) > pd_offset-0.060 & m.nodes(:,3) < pd_offset-0.035; % cut off lobes
    f_seplobes2_p(id_rimnodes)=1; % separate lobes from tube
    id_rimnodes = m.nodes(:,3) > pd_offset-0.085 & m.nodes(:,3) < pd_offset-0.065; % cut off lobes
    f_seplobes3_p(id_rimnodes)=1; % separate lobes from tube
    basenodes2 = m.nodes(:,3)< -0.17 & m.nodes(:,3) > -0.2;
    f_sepbase_p(basenodes2) = 1; % separate tube from base
    %m=leaf_set_seams(m,f_seam_p);
    % use f_seam as a marker
    
    
    
else
    if realtime>=240 && realtime<250 % over this initialisation phase
        % shrink the lobes (part of the startup process that ensures
        %triangles will not get too elongated during growth
        shrinkrate=-0.01;
        kapar_p=zeros(size(kapar_l));
        kbpar_p=zeros(size(kapar_l));
        kaper_p=zeros(size(kapar_l));
        kbper_p=zeros(size(kapar_l));
        kapar_p=shrinkrate.*pro(15,id_dtl_l.*inh(0.5, id_med_l));
        kbpar_p=shrinkrate.*pro(20,id_dtl_l.*inh (0.5, id_med_l));
    end
    disp(sprintf('\n cyc=%d dich=%d div=%d cenorg=%d v_dichdistorg=%d\n',...
        id_cyc_a,id_dich_a,id_div_a,id_cenorg_a,v_dichdistorg_a));
    %Induce sectors at 340hrs
    if (340>realtime-dt) && (340<realtime+dt)
        %         m = leaf_makesecondlayer( m, ...  % This function adds biological cells.
        %             'mode', 'each', ...  % Make biological cells randomly scattered over the flower.
        %             'relarea', 1/16000, ...   % Each cell has area was 1/16000 of the initial area of the flower.
        %             'probpervx', 'V_FLOWER', ... % induce transposed cells over whole corolla
        %             'numcells',4500,...%number of cells (that will become clones)
        %             'sides', 6, ...  % Each cell is approximated as a 6-sided regular polygon.
        %             'colors', [0.5 0.5 0.5], ...  % Default colour is gray but
        %             'allowoverlap', false, ...
        %             'colorvariation',1,... % Each cell is a random colour
        %             'add', true );  % These cells are added to any cells existing alread
    end
    
    
    % Setup should be common to all the models
    % Start setup
    if abs(realtime - 235) <0.5*dt %At 240 hrs
        %id_prox is fixed according to s_prox levels
        id_prox_p = 2*s_prox_l;
        id_prox_l = id_prox_p .* id_prox_a;
        
        %expression domain of id_lip established
        id_lip_p=zeros(size(id_lip_p));
        id_lipnodes = m.nodes(:,3) > pd_offset-0.095 & m.nodes(:,3) < pd_offset+0.01;
        id_lip_p(id_lipnodes) = 1;
        id_lip_l = id_lip_p * id_lip_a;
        m.morphogenclamp(id_lipnodes,id_lip_i) = 1;
        
        %expression domain of id_uptube established
        id_uptubenodes = m.nodes(:,3) > pd_offset-0.085 & m.nodes(:,3) < pd_offset-0.055;
        id_uptube_p(id_uptubenodes) = 1;
        id_uptube_l = id_uptube_p * id_uptube_a;
        
        %expression domain of id_lpb established
        id_lpbnodes = m.nodes(:,3) > pd_offset -0.035 & m.nodes(:,3) < pd_offset+0.01;
        id_lpb_p(id_lpbnodes) = 1;
        m.morphogenclamp(id_lpbnodes,id_lpb_i) = 1;
        
        %expression domain of id_mlobe established
        id_mlobenodes = m.nodes(:,3) > pd_offset -0.005 & m.nodes(:,3) < pd_offset+0.03;
        id_mlobe_p(id_mlobenodes) = 1;
        m.morphogenclamp(id_mlobenodes,id_mlobe_i) = 1;
        
        %expression domain of id_rim established
        id_rim_p=zeros(size(id_rim_p));
        id_rimnodes = m.nodes(:,3) > pd_offset-0.055 & m.nodes(:,3) < pd_offset-0.035;
        id_rim_p(id_rimnodes) = 1;
        id_rim_l = id_rim_p * id_rim_a;
        m.morphogenclamp(id_rimnodes,id_rim_i) = 1;
        
        %expression domain of id_tube established
        id_tube_p=zeros(size(id_tube_p));
        id_tubenodes = m.nodes(:,3)<= pd_offset-0.055 & m.nodes(:,3) > - 0.18;
        id_tube_p(id_tubenodes) = 1;
        id_tube_l = id_tube_p * id_tube_a;
        m.morphogenclamp(id_tubenodes,id_tube_i) = 1;
        
        %expression domain of id_lobe established
        id_lobe_p=zeros(size(id_lobe_p));
        id_lobenodes = m.nodes(:,3) >= pd_offset -0.055;%-0.035
        id_lobe_p(id_lobenodes) = 1;
        id_lobe_l = id_lobe_p * id_lobe_a;
        m.morphogenclamp(id_lobenodes,id_lobe_i) = 1;
        
        %XANA_expression domain of id_plt established
        id_plt_p=zeros(size(id_plt_p));
        id_pltnodes = m.nodes(:,3) > pd_offset-0.095 & m.nodes(:,3) < pd_offset-0.055; % it used to be 0.015 at upper boundary
        id_plt_p(id_pltnodes) = 1;
        id_plt_l = id_plt_p * id_plt_a;
        m.morphogenclamp(id_pltnodes,id_plt_i) = 1;
        
        %XANA_expression domain of id_lipcliff established
        id_lipcliff_p=zeros(size(id_lipcliff_p));
        %id_lipcliffnodes = m.nodes(:,3) > pd_offset-0.055 & m.nodes(:,3) < pd_offset-0.025; % it used to be 0.015 at upper boundary
        id_lipcliffnodes = m.nodes(:,3) > pd_offset-0.045 & m.nodes(:,3) < pd_offset-0.015; % it used to be 0.015 at upper boundary
        id_lipcliff_p(id_lipcliffnodes) = 1;
        id_lipcliff_l = id_lipcliff_p * id_lipcliff_a;
        m.morphogenclamp(id_lipcliffnodes,id_lipcliff_i) = 1;
        
        %XANA_expression domain of id_lipbend established
        id_lipdistal_p=zeros(size(id_lipdistal_p));
        %id_lipdistalnodes = m.nodes(:,3) > pd_offset -0.027 & m.nodes(:,3) < pd_offset+0.01; % it used to be -0.019 at lower boundary
        id_lipdistalnodes = m.nodes(:,3) > pd_offset -0.015 & m.nodes(:,3) < pd_offset+0.01; % it used to be -0.019 at lower boundary
        id_lipdistal_p(id_lipdistalnodes) = 1;
        id_lipdistal_l = id_lipdistal_p * id_lipdistal_a;
        m.morphogenclamp(id_lipdistalnodes,id_lipdistal_i) = 1;
        
    end
    % s_prox is activated by id_prox
    m.mgen_production(:,s_prox_i) = 0.1 * id_prox_l;
    
    % s_lpb is activated by id-lpb
    m.mgen_production(:,s_lpb_i) = 0.1 * id_lpb_l.*id_med_l;
    
    % s_lat is activated by id_lat
    m.mgen_production(:,s_lat_i) = 0.1 * id_lat_l;
    
    %s_med is activated by id_med
    m.mgen_production(:,s_med_i) = 0.1 * id_med_l; %
    
    
    if abs(realtime - 245) <0.5*dt % At 245 hrs
        %id_lat is fixed according to s_lat levels
        id_lat_p = 2.5*s_lat_l;
        
        %id_med is fixed according to s_med levels
        id_med_p = 2.5*s_med_l;
        %
        %             id_plt_p=id_tube_l.*id_lip_l;
        %             id_plt_l=id_plt_p.*id_plt_a;
        
        id_lip_p=id_lobe_l.*id_lip_l;
        id_lip_l=id_lip_p.*id_lip_a;
        
        m.userdata.extra_regions_ready=true;
        
        id_dtl_p=id_lobe_l.*inh(10,id_lip_l);
        id_dtl_l=id_dtl_p.*id_dtl_a;
        % set up the seams ready for dissection
        % (adult petals are dissected and flattened - Figures 5 and 6)
        id_rimnodes = m.nodes(:,3) > pd_offset-0.01 & m.nodes(:,3) < pd_offset+0.01; % cut off lobes
        
        %XANA_S_RIM
        s_rim_p(:) = 0;
        s_rim_p(:) = id_rim_p;
        m = leaf_mgen_conductivity( m, 'S_RIM', 0.00005 );
        m = leaf_mgen_absorption( m, 'S_RIM', 0.1);
        
        %XANA_S_SECVEIN
        s_secvein_p(:) = 0;
        s_secvein_p(:) = id_secvein_p;
        m = leaf_mgen_conductivity( m, 'S_SECVEIN', 0.0005 );
        m = leaf_mgen_absorption( m, 'S_SECVEIN', 0.1);
        
        f_seplobes1_p(id_rimnodes)=1; % separate lobes from tube
        %m=leaf_set_seams(m,f_seam_p);
        ind_mid=s_lat_p<0.18;
        f_seam2_p(ind_mid)=1; % all elements around mid
        ind_mid=s_lat_p<0.10;
        f_seam3_p(ind_mid)=1; % all elements around mid
        ind_vent_lat=s_rad_p>0.02 & s_rad_p<0.06;
        f_seam4_p(ind_vent_lat)=1; % all elements ventral side of laterals
        ind_vent_lat=s_rad_p<=0.04;
        f_seamventral_p(ind_vent_lat)=1;
        ind_dorsal_lat=s_rad_p>0.47;
        f_seamdorsal_p(ind_dorsal_lat)=1;
    end
    
    if abs(realtime - 250) <0.5*dt % At 245 hrs
        m = leaf_mgen_conductivity( m, 'S_RIM', 0 );
        m = leaf_mgen_absorption( m, 'S_RIM', 0 );
        m = leaf_mgen_conductivity( m, 'S_SECVEIN', 0 );
        m = leaf_mgen_absorption( m, 'S_SECVEIN', 0);
        
        id_flank_p(:) = 10*(s_secvein_p .* s_rim_p .*s_lat_p);
        
    end
    
    v_dichdistorg_p(:)=1; % visualisation
    % End setup
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%  START INDIVIDUAL MODELS %%%%%%%%%%%%%%%%%%%%%%
    
    if strcmpi(modelnumber,'Fig 9B Model 7')
        % P EQUATIONS
        % START POLARISER section (it is slightly mixed with gene networks)
        m = leaf_mgen_conductivity( m, 's_distorg', 0.002); %was 0.1 distorg_diff);% diffusion constant
        m = leaf_mgen_absorption( m, 's_distorg', 0.05);  % it will not decay everywhere
        m.globalProps.mingradient=0.01; % this ensures that all polariser gradients are frozen
        m = leaf_mgen_conductivity( m, 'Polariser', 0.002);% diffusion constant WAS 0.001
        m = leaf_mgen_dilution( m, 'Polariser', false );% it will not dilute with growth
        m = leaf_mgen_absorption( m, 'Polariser', 0.1);  % it will not decay everywhere
        
        m = leaf_mgen_conductivity( m, 's_cenorg', 0.002); %was 0.1 distorg_diff);% diffusion constant
        m = leaf_mgen_absorption( m, 's_cenorg', 0.05);  % it will not decay everywhere
        
        %\latex \subsection{PRN (changes)}
        % POLARISER PRODUCTION
        if realtime>230 % From time 230  generate polariser
            id_proxorg_p = (id_prox_l >0.8);                                             % Eqn 2  CHECKED
            id_proxorg_l = id_proxorg_p * id_proxorg_a;
        end
        
        % POLARISER DECAY
        % switch on distorg with id_dist in the absence of id_lat
        if abs(realtime - 230) <0.5*dt % At time 230 hrs
            %                 id_distorg_p = 0.25* id_dist_l .* inh (100, id_lat_l);%  CHECKED % Eqn 3 (used in production equation)
            %                 id_distorg_l = id_distorg_p * id_distorg_a;
            %                 m.userdata.all_regions_ready=true;
            id_distorg_p = 0.25* id_dist_l .* inh (100, id_lat_l)+ 0.1*id_dich_l .* id_dist_l .* id_cyc_l;           % Eqn 3 (used in production equation)
            id_distorg_l = id_distorg_p * id_distorg_a;
            m.userdata.all_regions_ready=true;
            v_dichdistorg_p=v_dichdistorg_p.*id_dich_l.* id_dist_l .* id_cyc_l;
        end
        
        if abs(realtime - 245) <0.5*dt %At 240 hrs
            %id_cenorg is activated by id_rim, id_lat and
            % id_div in absence of id_rad
            id_cenorg_p = 0.6*id_lat_l .* id_rim_l .* id_div_l ;                                   % Eqn 4 % CHECKED
            id_cenorg_l = id_cenorg_p * id_cenorg_a;
        end
        
        m.mgen_production(:,polariser_i) = 0.1 * (id_proxorg_l + 1 - P .* (id_distorg_l+3*id_late_p.*id_cenorg_l)); % Eqn 19 decay part of WAS S_DISTORG BY MISTAKE
        %m.mgen_production(:,polariser_i) = m.mgen_production(:,polariser_i) - 0.1 * P .* (id_distorg_l+3*id_late_p.*id_cenorg_l);  % Eqn 19 decay part of WAS S_DISTORG BY MISTAKE
        
        % LATER MODULATION OF DISTORG AND DIST
        % s_dist is activated by id_dist
        m.mgen_production(:,s_dist_i) = 0.1 * id_dist_l.* inh (100, id_lat_l);          % Eqn 17  CHECKED
        
        % s_distorg is activated by id_distorg
        m.mgen_production(:,s_distorg_i) = 0.1 * id_distorg_l ;                         % Eqn 15 (but  CHECKED
        
        % s_cenorg is activated by id_cenorg
        m.mgen_production(:,s_cenorg_i) = 0.1 * id_cenorg_l;                             % Eqn 16 % CHECKED
        
        if abs(realtime - 340) <0.5*dt % At time 340 hrs
            id_cenorg_p = s_cenorg_l;                                                   % Eqn 11 % CHECKED
            id_cenorg_l = id_cenorg_p .* id_cenorg_a;
            id_distorg_p = s_distorg_l;                                                 % Eqn 12  CHECKED
            id_distorg_l = id_distorg_p .* id_distorg_a;
            id_dist_p = s_dist_l;                                                       % Eqn 13  CHECKED (NOT NEEDED IN THIS MODEL)
            id_dist_l = id_dist_p .* id_dist_a;
        end
        % End Polariser
        %\latex \subsection{GRN (changes)}
        
        % id_rad is activated by id_cyc and id_dich at ALL TIMES
        id_rad_p = id_cyc_l +  0.3* id_dich_l;                                           % Eqn 1 % CHECKED
        id_rad_l = id_rad_p * id_rad_a;
        
        % s_rad is activated by id_rad
        m.mgen_production(:,s_rad_i) = 0.1 * id_rad_l .* pro (0.5, id_tube_l);          % Eqn 14 % CHECKED
        
        % I suspect that it is the order within this set of equations that matters because
        % cenorg depends on div and rim in the setup secton that also happen
        % at 240
        if abs(realtime - 240) <0.5*dt %At 240 hrs
            % continue with the initial regionalisation
            % replace [0,1] values with reals for quantitative growth
            
            %switch on id_div
            id_div_p=inh(100, id_rad_l);                           % Eqn 5 % CHECKED
            id_div_l=id_div_p.*id_div_a;
            
            %id_laterals is activated by id_div
            id_lts_p(:) = id_div_l;                                                      % Eqn 6 % CHECKED
            id_lts_l=id_lts_p * id_lts_a;
        end
        
        if abs(realtime - 250) <0.5*dt % At time 250 hrs
            %id_cyc and id_dich are inhibited by id_div in absence of id_rad.
            %This is to ensure that in a id_rad mutant,
            %id_cyc and id_dich are inactivated. (Also see constitutive cyc
            %below)
            id_cyc_p = id_cyc_l.* inh (4, id_div_l );                                    % Eqn 7 % CHECKED
            id_cyc_l = id_cyc_p * id_cyc_a;
            
            id_dich_p = id_dich_l .* inh (4, id_div_l );                                 % Eqn 8 % CHECKED
            id_dich_l = id_dich_p * id_dich_a;
            
            % one shot setup regionalising lateral petals
            id_lts_p=5*id_lts_l.* s_rad_l;                                               % Eqn 9 % CHECKED
            id_lts_l=id_lts_p.*id_lts_a;
        end
        if abs(realtime - 340) <0.5*dt % At time 340 hrs
            % narrow domain of div expression with id_laterals and id_rad at ~ 14 days
            id_div_p = id_div_l .* inh (5, id_lts_l);                                   % Eqn 10
            id_div_l = id_div_p.*id_div_a;
        end
        
        %             if realtime >= 352.5
        %                 m = leaf_setpolfrozen (m, id_late_p, id_late_p);
        % %                 m = leaf_setpolfrozen (m, id_lip_p, id_lip_p);
        % %                 m = leaf_setpolfrozen (m, id_tube_p, id_tube_p);
        %             end
        
        % End Gene Networks
        
        %\latex \subsection{KRN (changes)}
        if realtime >= 250 && realtime<570  % Between time 250 570
            % Compute growth rates
            
            % Note: the factors are vectors, one element per vertex (node)
            % as a result, with each statement, different things happen in
            % different regions of the model.
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            Kpar=... % Eqn 260
                0.013 ... % Background growth parallel to the polariser  Eqn 200
                .* inh (0.2, id_prox_l)... %keeps base small                              Eqn 24 % CHECKED
                .* pro (0.4,id_tube_l .* inh (100, id_plt_l))... %promote growth of tube  Eqn 25 % CHECKED
                .* pro (1.4, id_early_p .* id_lip_l .*inh (100, id_rad_l.* inh(100, id_dich_l.* id_lat_l)))... %form lower face   Eqn 30 % CHECKED
                .* inh (10,  id_early_p .* id_rim_l )... %inhibit rim to help with bending back    Eqn 32 % CHECKED
                .* inh (0.5, id_late_p .* id_med_l.* id_lobe_l)... %reduce growth of midlobe   Eqn 38  % CHECKED
                .* pro (2.4, id_early_p .* id_div_l.*id_plt_l )... %make lower tube arch over      Eqn 31 % CHECKED
                .* inh (1, id_rad_l .*id_lip_l .* inh(3, id_dich_l .* id_lat_l))... %prevent lip growing at edge of dorsals   Eqn 26 % CHECKED
                .* inh (1, id_rad_l .*id_plt_l .* inh(10, id_dich_l) .* inh (20, id_cyc_l .* inh (40, id_lat_l.^2)))... %prevent palate growing at edge of dorsals   Eqn 27 % CHECKED
                .* pro (0.3, id_cyc_l .* id_dtl_l  .* inh (0.5, id_dich_l))... %promote growth of lateral dorsal lobe         Eqn 28 % CHECKED
                .* pro (0.45, (id_cyc_l + 0.2* id_dich_l) .* id_plt_l)... %promote medial palate growth with CYC              Eqn 29 % CHECKED
                .* pro (2.2, id_late_p .* id_lts_l .* id_med_l .* inh (0.5, id_lpb_l).* inh(4, id_lat_l) .*(id_lip_l + 0.3* id_plt_l))...%  promote lateral growth   Eqn 42 % CHECKED
                ... %%%%%%% FINE TUNING %%%%%% better formed palate
                .* inh (10,  id_late_p .* id_div_l .*id_cenorg_l.* inh(2, id_lpb_l).* pro (10, id_plt_l .* id_med_l))...% prevent lips protruding   Eqn 39
                .* pro (0.2, id_late_p .* id_div_l .* id_plt_l)... %promote growth of palate   Eqn 40
                .* inh (0.7, id_late_p .* id_lts_l .* id_lat_l .* (id_lip_l + id_plt_l)); %inhibit growth at edge of laterals   Eqn 41
            
            
            
            %                 DKapar=1 ... %%%%%%% FINE TUNING %%%%%%
            %                     .* inh (1, id_dich_l.* id_lat_l  .* id_lpb_l)...                               % Eqn 33
            %                     .* inh (1, id_early_p .* (id_cyc_l + id_div_l) .*id_mlobe_l .* inh (5, id_dich_l))... %bend up distal lobes   Eqn 34
            %                     .* pro (1, id_early_p .*  id_div_l .* id_lpb_l.* inh (5, id_med_l))... % bend out lips   Eqn 35
            
            %AKapar=0.05* id_early_p .* id_rim_l .* ( 1.2*id_div_l + 0.3* id_lts_l.* id_med_l  + 0.5*(id_cyc_l +id_dich_l)); % bend back petals   Eqn 36 % CHECKED
            
            %%%%%%% FINE TUNING %%%%%%
            %                 DKbpar=pro (2, id_dich_l .* id_lat_l  .* id_lpb_l); % bend dorsal lip forward   Eqn 37  % CHECKED
            
            % Background growth normal to the polariser   Eqn 201
            
            Kper=... % Eqn 250
                0.0075 ... % Background growth normal to the polariser   Eqn 201
                .* inh (0.2, id_prox_l)... %keep base small   Eqn 43   % CHECKED
                .* pro(2, id_early_p .* id_dist_l .* inh (20, id_lat_l))... %  widen distal lobe    Eqn 46  % CHECKED
                .* inh (0.3, id_late_p .* id_lobe_l .* id_med_l)... %keep medial lobe narrow   Eqn 50  % CHECKED
                .* pro (1,   id_late_p .* id_dist_l .* pro (1.2, (id_cyc_l + id_div_l)))... % increase width of distal lobe   Eqn 49  % CHECKED
                .* inh (1.3, id_late_p .* id_div_l .* inh(2.5, id_lobe_l.* inh (2,id_lpb_l)).* pro(1, id_plt_l))...  %narrow ventral petals    Eqn 51 %%%%%%% FINE TUNING %%%%%%
                .* pro (0.1, id_cyc_l  .* pro(1.5,id_lip_l))... % increase width of dorsal petal    Eqn 44  % CHECKED
                .* inh (6, id_early_p .* id_lts_l.*id_med_l )... %  keep laterals narrow initially   Eqn 45
                .* pro (0.2, id_late_p .* id_lts_l .* id_med_l .*  id_plt_l); %promote growth of lateral width   Eqn 52   ADDD Fig 6G Model 3
            
            %%%%%%% FINE TUNING %%%%%%
            %                 DKaper = pro (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l));%bend out lips;   Eqn 47
            %                 DKbper = inh (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l));%bend out lips;   Eqn 48
            
            
            % compute the growth rates
            kapar=Kpar ...% Eqn 300
                .* inh (1, id_dich_l.* id_lat_l  .* id_lpb_l)...                               % Eqn 33
                .* inh (1, id_early_p .* (id_cyc_l + id_div_l) .*id_mlobe_l .* inh (5, id_dich_l))... %bend up distal lobes   Eqn 34
                .* pro (1, id_early_p .*  id_div_l .* id_lpb_l.* inh (5, id_med_l))... % bend out lips   Eqn 35
                +0.05* id_early_p .* id_rim_l .* ( 1.2*id_div_l + 0.3* id_lts_l.* id_med_l  + 0.5*(id_cyc_l +id_dich_l)); % Eqn 301
            kbpar=Kpar ... % Eqn 310
                .*pro (2, id_dich_l .* id_lat_l  .* id_lpb_l); % Eqn 37
            kaper=Kper ... % Eqn 320
                .* pro (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l));%bend out lips;   Eqn 47
            kbper=Kper ... % Eqn 330
                .*inh (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l));%bend out lips;   Eqn 48
            knor =0.003; % thickness tracks model using Scaled USE KNOR(anti-curl) Eqn 305
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif strcmpi(modelnumber,'Figure 9_Supplement 1A. Ground state')
        % P EQUATIONS
        % START POLARISER section (it is slightly mixed with gene networks)
        m = leaf_mgen_conductivity( m, 's_distorg', 0.002); %was 0.1 distorg_diff);% diffusion constant
        m = leaf_mgen_absorption( m, 's_distorg', 0.05);  % it will not decay everywhere
        m.globalProps.mingradient=0.01; % this ensures that all polariser gradients are frozen
        m = leaf_mgen_conductivity( m, 'Polariser', 0.002);% diffusion constant WAS 0.001
        m = leaf_mgen_dilution( m, 'Polariser', false );% it will not dilute with growth
        m = leaf_mgen_absorption( m, 'Polariser', 0.1);  % it will not decay everywhere
        
        m = leaf_setproperty (m, ...
            'twosidedpolarisation', true, ... % the mesh has two-sided polarisation
            'mingradient', 0.0, ... % No threshold for freezing the polariser
            'usedpolfreezebc', true); % a rather subtle choice of a feature of the behaviour of frozen gradients. It probably doesn't make much difference.
        
        m = leaf_setpolfrozen (m, false); %Initially the gradient is not frozen anywhere
        
        %\latex \subsection{PRN (changes)}
        % POLARISER PRODUCTION
        if realtime>230 % From time 230  generate polariser
            id_proxorg_p = (id_prox_l >0.8);                                             % Eqn 2  CHECKED
            id_proxorg_l = id_proxorg_p * id_proxorg_a;
        end
        
        % POLARISER DECAY
        % switch on distorg with id_dist in the absence of id_lat
        if abs(realtime - 230) <0.5*dt % At time 230 hrs
            id_distorg_p = 0.25* id_dist_l .* inh (100, id_lat_l)+ 0.1*id_dich_l .* id_dist_l .* id_cyc_l;           % Eqn 3 (used in production equation)
            id_distorg_l = id_distorg_p * id_distorg_a;
            m.userdata.all_regions_ready=true;
            v_dichdistorg_p=v_dichdistorg_p.*id_dich_l.* id_dist_l .* id_cyc_l;
        end
        
        if abs(realtime - 240) <0.5*dt %At 240 hrs
            
            %S_RAD NEW Rebocho et al%
            s_rad_p(:) = 0;
            s_rad_p(:) = 2*id_rad_p (:);
            m = leaf_mgen_conductivity( m, 'S_RAD', 0.002 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0.1 );
        end
        
        m.mgen_production(:,polariser_i) = 0.1 * (id_proxorg_l + 1 - P .* id_distorg_l); % Eqn 19 decay part of WAS S_DISTORG BY MISTAKE
        
        % LATER MODULATION OF DISTORG AND DIST
        % s_dist is activated by id_dist
        m.mgen_production(:,s_dist_i) = 0.1 * id_dist_l.* inh (100, id_lat_l);          % Eqn 17  CHECKED
        
        % s_distorg is activated by id_distorg
        m.mgen_production(:,s_distorg_i) = 0.1 * id_distorg_l ;                         % Eqn 15 (but  CHECKED
        
        if abs(realtime - 340) <0.5*dt % At time 340 hrs
            id_distorg_p = s_distorg_l;                                                 % Eqn 12  CHECKED
            id_distorg_l = id_distorg_p .* id_distorg_a;
            id_dist_p = s_dist_l;                                                       % Eqn 13  CHECKED (NOT NEEDED IN THIS MODEL)
            id_dist_l = id_dist_p .* id_dist_a;
        end
        % End Polariser
        %\latex \subsection{GRN (changes)}
        
        % id_rad is activated by id_cyc and id_dich at ALL TIMES
        id_rad_p = id_cyc_l +  0.3* id_dich_l;                                           % Eqn 1 % CHECKED
        id_rad_l = id_rad_p * id_rad_a;
        
        if abs(realtime - 250) <0.5*dt % At time 250 hrs
            %id_cyc and id_dich are inhibited by id_div in absence of id_rad.
            %This is to ensure that in a id_rad mutant,
            %id_cyc and id_dich are inactivated. (Also see constitutive cyc
            %below)
            id_cyc_p = id_cyc_l.* inh (4, id_div_l );                                    % Eqn 7 % CHECKED
            id_cyc_l = id_cyc_p * id_cyc_a;
            
            id_dich_p = id_dich_l .* inh (4, id_div_l );                                 % Eqn 8 % CHECKED
            id_dich_l = id_dich_p * id_dich_a;
            
            %S_RAD gradient stabilises- NEW Rebocho et al%
            m = leaf_mgen_conductivity( m, 'S_RAD', 0 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0 );
            
            %switch on id_div_NEW Rebocho et al%
            id_div_p=zeros(size(id_div_p));
            id_div_p=inh(100, id_rad_l).* inh (10, s_rad_p);                           % Eqn 5 % CHECKED
            id_div_l=id_div_p.*id_div_a;
            
            id_lipvisual_p = id_lip_l.*inh(100, id_rim_l);
            
        end
        
        if realtime >= 250 && realtime<352.5 %Freeze polarity in the tube region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, id_tube_p, id_tube_p);
        end
        
        if realtime >= 352.5 %Freeze polarity in the tube and Lip region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, (id_lip_p + id_tube_l), (id_lip_p + id_tube_l));
        end
        
        % End Gene Networks
        
        %\latex \subsection{KRN (changes)}
        if realtime >= 250 && realtime<570  % Between time 250 570
            % Compute growth rates
            
            % Note: the factors are vectors, one element per vertex (node)
            % as a result, with each statement, different things happen in
            % different regions of the model.
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % need to delete effect of cenorg on polarity
            
            Kpar=... % Eqn 260
                0.013 ... % Background growth parallel to the polariser  Eqn 200
                ...%%%%%equations maintained from snapdragon published Fig 9B model7%%%%%%
                .* inh (0.2, id_prox_l)... %keeps base small                              Eqn 24
                .* pro (0.4,id_tube_l .* inh (100, id_plt_l))... %promote growth of tube  Eqn 25
                .* pro (1.4, id_early_p .* id_lip_l .*inh (100, id_rad_l.* inh(100, id_dich_l.* id_lat_l)))...  To grow the dorsal lip a bit more   Eqn 30
                .* inh (10,  id_early_p .* id_rim_l )... %inhibit rim to help with bending back    Eqn 32
                .* inh (0.5, id_late_p .* id_med_l.* id_lobe_l)... %  Eqn 38
                .* inh (1, id_rad_l .*id_lip_l .* inh(6, id_dich_l .* id_lat_l))... % prevent lip growing at edge of dorsals  Eqn 26
                .* inh (1, id_rad_l .*id_plt_l .* inh (15, id_dich_l .* inh (5, id_lat_l) ) .* inh (30, id_cyc_l .* inh (40, id_lat_l.^2)))... %prevent palate growing at edge of dorsals   Eqn 27
                .* pro (0.3, id_cyc_l .* id_dtl_l  .* inh (0.5, id_dich_l))... %promote growth of lateral dorsal lobe         Eqn 28
                .* pro (0.45, (id_cyc_l + 0.2* id_dich_l) .* id_plt_l)... %promote medial palate growth with CYC              Eqn 29
                ...%%%%%%equations removed from published Fig 9B model7
                ...%.* pro (2.4, id_early_p .* id_div_l.*id_plt_l )... %make lower tube arch over      Eqn 31
                ...% .* pro (2.2, id_late_p .* id_lts_l .* id_med_l .* inh (0.5, id_lpb_l).* inh(4, id_lat_l) .*(id_lip_l + 0* id_plt_l));%  promote lateral growth   Eqn 42
                
            Kper=... % Eqn 250
                0.0075 ... % Background growth normal to the polariser   Emqn 201
                ...%%%%%equations maintained from snapdragon published Fig 9B odel7%%%%%%
                .* inh (0.2, id_prox_l)... %keep base small   Eqn 43
                .* pro (2, id_early_p .* id_dist_l .* inh (20, id_lat_l))... %  widen distal lobe    Eqn 46
                .* inh (0.3, id_late_p .* id_lobe_l .* id_med_l)... %keep medial lobe narrow   Eqn 50
                .* pro (1, id_late_p .* id_dist_l .* pro (1.2, (id_cyc_l + 0*id_div_l)))... %  increase width of distal lobe  Eqn 49
                .* pro (0.1, id_cyc_l  .* pro(1.5,id_lip_l)); % increase width of dorsal petal    Eqn 44
            ...%%%%%%equations removed from published Fig 9B model7
                ....* inh (1.3, id_late_p .* id_div_l .* inh(2.5, id_lobe_l.* inh (2,id_lpb_l)).* pro(1, id_plt_l))...  %narrow ventral petals    Eqn 51 %%%%%%% FINE TUNING %%%%%%
                ...%.* inh(6, id_early_p .* id_lts_l.*id_med_l )... %  keep laterals narrow initially   Eqn 45
                ...%.* pro (0.2, id_late_p .* id_lts_l .* id_med_l .*  id_plt_l); %promote growth of lateral width   Eqn 52   ADDD Fig 6G Model 3
                
            % compute the growth rates
            kapar=Kpar ...% Eqn 300
                .* inh (1, id_dich_l.* id_lat_l  .* id_lpb_l)...  % bend dorsal petals Eqn 33
                .* inh (1, id_early_p .*(id_cyc_l + 0*id_div_l) .*id_mlobe_l .* inh (5, id_dich_l))...  %bend up distal lobes   Eqn 34
                ....* inh (0.2, id_div_l .* id_lpb_l.* inh (5, id_med_l)) ... %  bend out lips Eqn 35
                +0.05* id_early_p .* id_rim_l .* ( 0*id_div_l + 0* id_lts_l.* id_med_l  + 0.5*(id_cyc_l +id_dich_l)); % Eqn 301
            kbpar=Kpar... % Eqn 310
                .* pro (2, id_dich_l .* id_lat_l  .* id_lpb_l); %Eqn 37 promote abaxial growth of lip triangle in between dorsal petal
            kaper=Kper; % Eqn 320
            ....* pro (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l)); %bend out lips Eqn 47
                kbper=Kper; % Eqn 330
            ....* inh (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l));%bend out lips  Eqn 48
                knor =0.003;
            v_kaniso_p = log( kbpar_p./kbper_p);% visualise growth anisotropy in ln scale
            v_karea_p = kbper_p + kbpar_p;% visualise areal growth rate
            
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    elseif strcmpi(modelnumber,'Figure 9K.  div mutant')
       % P EQUATIONS
        % START POLARISER section (it is slightly mixed with gene networks)
        m = leaf_mgen_conductivity( m, 's_distorg', 0.002); %was 0.1 distorg_diff);% diffusion constant
        m = leaf_mgen_absorption( m, 's_distorg', 0.05);  % it will not decay everywhere
        m.globalProps.mingradient=0.01; % this ensures that all polariser gradients are frozen
        m = leaf_mgen_conductivity( m, 'Polariser', 0.002);% diffusion constant WAS 0.001
        m = leaf_mgen_dilution( m, 'Polariser', false );% it will not dilute with growth
        m = leaf_mgen_absorption( m, 'Polariser', 0.1);  % it will not decay everywhere
        
        m = leaf_setproperty (m, ...
            'twosidedpolarisation', true, ... % the mesh has two-sided polarisation
            'mingradient', 0.0, ... % No threshold for freezing the polariser
            'usedpolfreezebc', true); % a rather subtle choice of a feature of the behaviour of frozen gradients. It probably doesn't make much difference.
        
        m = leaf_setpolfrozen (m, false); %Initially the gradient is not frozen anywhere
        
        %\latex \subsection{PRN (changes)}
        % POLARISER PRODUCTION
        if realtime>230 % From time 230  generate polariser
            id_proxorg_p = (id_prox_l >0.8);                                             % Eqn 2  CHECKED
            id_proxorg_l = id_proxorg_p * id_proxorg_a;
        end
        
        % POLARISER DECAY
        % switch on distorg with id_dist in the absence of id_lat
        if abs(realtime - 230) <0.5*dt % At time 230 hrs
            id_distorg_p = 0.25* id_dist_l .* inh (100, id_lat_l)+ 0.1*id_dich_l .* id_dist_l .* id_cyc_l;           % Eqn 3 (used in production equation)
            id_distorg_l = id_distorg_p * id_distorg_a;
            m.userdata.all_regions_ready=true;
            v_dichdistorg_p=v_dichdistorg_p.*id_dich_l.* id_dist_l .* id_cyc_l;
        end
        
        if abs(realtime - 240) <0.5*dt %At 240 hrs
            
            %S_RAD NEW Rebocho et al%
            s_rad_p(:) = 0;
            s_rad_p(:) = 2*id_rad_p (:);
            m = leaf_mgen_conductivity( m, 'S_RAD', 0.002 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0.1 );
        end
        
        m.mgen_production(:,polariser_i) = 0.1 * (id_proxorg_l + 1 - P .* id_distorg_l); % Eqn 19 decay part of WAS S_DISTORG BY MISTAKE
        
        % LATER MODULATION OF DISTORG AND DIST
        % s_dist is activated by id_dist
        m.mgen_production(:,s_dist_i) = 0.1 * id_dist_l.* inh (100, id_lat_l);          % Eqn 17  CHECKED
        
        % s_distorg is activated by id_distorg
        m.mgen_production(:,s_distorg_i) = 0.1 * id_distorg_l ;                         % Eqn 15 (but  CHECKED
        
        if abs(realtime - 340) <0.5*dt % At time 340 hrs
            id_distorg_p = s_distorg_l;                                                 % Eqn 12  CHECKED
            id_distorg_l = id_distorg_p .* id_distorg_a;
            id_dist_p = s_dist_l;                                                       % Eqn 13  CHECKED (NOT NEEDED IN THIS MODEL)
            id_dist_l = id_dist_p .* id_dist_a;
        end
        % End Polariser
        %\latex \subsection{GRN (changes)}
        
        % id_rad is activated by id_cyc and id_dich at ALL TIMES
        id_rad_p = id_cyc_l +  0.3* id_dich_l;                                           % Eqn 1 % CHECKED
        id_rad_l = id_rad_p * id_rad_a;
        
        if abs(realtime - 250) <0.5*dt % At time 250 hrs
            %id_cyc and id_dich are inhibited by id_div in absence of id_rad.
            %This is to ensure that in a id_rad mutant,
            %id_cyc and id_dich are inactivated. (Also see constitutive cyc
            %below)
            id_cyc_p = id_cyc_l.* inh (4, id_div_l );                                    % Eqn 7 % CHECKED
            id_cyc_l = id_cyc_p * id_cyc_a;
            
            id_dich_p = id_dich_l .* inh (4, id_div_l );                                 % Eqn 8 % CHECKED
            id_dich_l = id_dich_p * id_dich_a;
            
            %S_RAD gradient stabilises- NEW Rebocho et al%
            m = leaf_mgen_conductivity( m, 'S_RAD', 0 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0 );
            
            %switch on id_div_NEW Rebocho et al%
            id_div_p=zeros(size(id_div_p));
            id_div_p=inh(100, id_rad_l).* inh (10, s_rad_p);                           % Eqn 5 % CHECKED
            id_div_l=id_div_p.*id_div_a;
            
            %Specification of hinge region_NEW Rebocho et al%
            id_hinge_p = (s_rad_l<0.39).* (s_rad_l>0.03);
            id_hinge_l = id_hinge_p .* id_hinge_a;
            
            %Specification of LP_DP region_NEW Rebocho et al%
            id_lpdp_p = (s_rad_l<0.55).* (s_rad_l>0.2).* id_lat_l;
            id_lpdp_l = id_lpdp_p .* id_lpdp_a;
            
            %id_msinus_NEW Rebocho et al%
            id_msinus_p = (id_lat_l>0.65).* id_lobe_l.* inh (10, s_rim_l) ; %graded expression of ID MSINUS from sinus up to rim
            id_msinus_l = id_msinus_p .* id_msinus_a;
            
            id_lipvisual_p = id_lip_l.*inh(100, id_rim_l);
            
            
        end
        
        if realtime >= 250 && realtime<352.5 %Freeze polarity in the tube region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, id_tube_p, id_tube_p);
        end
        
        if realtime >= 340 % Reorientation of polarity at the ventral-lateral sinuses_NEW Rebocho et al%
            P(id_msinus_p.* id_div_l >0.7) = 0.9;
            m.morphogenclamp((id_msinus_p.* id_div_l> 0.7), polariser_i) = 1;
            m = leaf_setpolfrozen (m, false, id_lip_p);
        end
        
        if realtime >= 352.5 %Freeze polarity in the tube and Lip region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, (id_lip_p + id_tube_l), (id_lip_p + id_tube_l));
        end
        
        % End Gene Networks
        
        %\latex \subsection{KRN (changes)}
        if realtime >= 250 && realtime<570  % Between time 250 570
            % Compute growth rates
            
            % Note: the factors are vectors, one element per vertex (node)
            % as a result, with each statement, different things happen in
            % different regions of the model.
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % need to delete effect of cenorg on polarity
           

Kpar=... % Eqn 260
                0.013 ... % Background growth parallel to the polariser  Eqn 200
                ...%%%%%equations maintained from snapdragon published Fig 9B model7%%%%%%
                .* inh (0.2, id_prox_l)... %keeps base small                              Eqn 24
                .* pro (0.4,id_tube_l .* inh (100, id_plt_l))... %promote growth of tube  Eqn 25
                .* pro (1.4, id_early_p .* id_lip_l .* id_dich_l.* id_lat_l)... %Modified for Rebocho et al - To grow the dorsal lip a bit more - modified from .* pro (1.4, id_early_p .* id_lip_l .*inh (100, id_rad_l.* inh(100, id_dich_l.* id_lat_l)))  Eqn 30
                .* inh (0.5, id_late_p .* id_med_l.* id_lobe_l.* id_cyc_l)... %Modified for Rebocho et al - to separate the control of the growth of dorsal lobe from the lower petal lobes - MODIFIED from .* inh (0.5, id_late_p .* id_med_l.* id_lobe_l) Eqn 38
                .* inh (0.3, id_late_p .* id_lobe_l.*pro (0.5, id_med_l).* inh (5, id_dist_l ).* inh (5,id_lipcliff_l).* inh (100, id_cyc_l))... %NEW Rebocho et al - TO GROW VENTRAL LOBE MORE
                .* inh (1, id_rad_l .*id_lip_l .* inh(6, id_dich_l .* id_lat_l))... %Modified for Rebocho et al - prevent lip growing at edge of dorsals  MODIFIED from .* inh (1, id_rad_l .*id_lip_l .* inh(3, id_dich_l .* id_lat_l) )Eqn 26
                .* inh (1, id_rad_l .*id_plt_l .* inh (15, id_dich_l .* inh (5, id_lat_l) ) .* inh (30, id_cyc_l .* inh (40, id_lat_l.^2)))... %prevent palate growing at edge of dorsals   Eqn 27 % CHECKED
                .* pro (0.3, id_cyc_l .* id_dtl_l  .* inh (0.5, id_dich_l))... %promote growth of lateral dorsal lobe         Eqn 28 % CHECKED
                .* pro (0.45, (id_cyc_l + 0.2* id_dich_l) .* id_plt_l)... %promote medial palate growth with CYC              Eqn 29 % CHECKED
                ...%%%%%% New added equations for Rebocho et al, for div mutant%%%%%
                .* pro (0.3, id_early_p .* id_lobe_l.* inh (100, id_lip_l).*inh (1, id_rad_l))...% Grow lobe region at early stages - necessary to ovatined bigger lobes in all petals
                .* pro (2, id_early_p .* id_lipdistal_l.*inh (5, s_rad_l) .*inh (100, id_rad_l).* inh (0.7, id_div_l))...% grow lipdistal only at ealy stages
                .* pro (0.5, id_late_p .* id_lipdistal_l .*inh (5, s_rad_l).*inh (100, id_rad_l).* inh (0.7, id_div_l))...% grow lipdistal only at late stages
                .* inh (2, id_late_p .*(s_rim_l.^2).*inh (100, id_rad_l))... %Keep rim region short
                .* pro (0.35, id_late_p .* id_dtl_p .* id_hinge_l .*inh (100, id_lip_l))... %grow distal lobes in hinge region
                .* inh (2, s_rad_l .*(id_lip_p + 0.5*id_uptube_l).* inh (100,s_rad_l > 0.25))...%Keep hinge region small
                ....%%%%%%New added equations for Rebocho et al, 2016 for wild-type%%%%%
                .* pro (2.5, id_early_p.* id_div_l.*id_lipcliff_p .*inh (100, id_rim_l))...% grow lipcliff only at early stages
                .* pro (0.2, id_late_p.* id_div_l.*id_lipcliff_p .*inh (100, id_rim_l))....% grow lipcliff only at late stages ...% GROW PLT EARLY IS ESSENTIAL FOR BENDING THEPETALS AND HELPING WITH SHAPE
                .* pro (2, id_early_p .* id_div_l .* id_plt_l )... .* pro (0.4, id_uptube_p))... % promote growth in palate with DIV
                .* pro (0.45, id_late_p .* id_div_l .* id_plt_l .* inh (2, id_uptube_l .* inh (100, id_lpdp_l)))... % promote growth in palate with DIV
                .* pro (0.8, id_late_p.* id_lat_l .* (id_lip_l + id_uptube_l ) .*inh (100, s_rad_l>0.4).* pro(1, id_lpdp_l).* inh (100, id_lip_l.*id_lpdp_l))...%vertical arms of orthogonal anisotropy
                .* inh (0.5, id_late_p.* id_med_l .* (id_lip_l + id_uptube_l) .* pro (2, s_rim_l) .*inh (100, id_rad_l))... % keep midvein region narrow
                .* inh (0.8, id_late_p.* id_med_l .* (0.5*id_plt_p + id_lipcliff_l).*(id_div_l>0.97)); %make ventral cleft
            ...%%%%%%equations removed from published Fig 9B model7
                ...%.* inh (10,  id_early_p .* id_rim_l )... %inhibit rim to help with bending back    Eqn 32
                ...%.* pro (2.4, id_early_p .* id_div_l.*id_plt_l )... %make lower tube arch over      Eqn 31
                ...% .* pro (2.2, id_late_p .* id_lts_l .* id_med_l .* inh (0.5, id_lpb_l).* inh(4, id_lat_l) .*(id_lip_l + 0* id_plt_l));%  promote lateral growth   Eqn 42
                
            Kper=... % Eqn 250
                0.0075 ... % Background growth normal to the polariser   Emqn 201
                ...%%%%%equations maintained from snapdragon published Fig 9B odel7%%%%%%
                .* inh (0.2, id_prox_l)... %keep base small   Eqn 43
                .* pro (2, id_early_p .* id_dist_l .* inh (20, id_lat_l))... %  widen distal lobe    Eqn 46
                .* pro (1, id_late_p .* id_dist_l .* pro (1.8, (id_cyc_l + id_div_l)))... % Modified for Rebocho et al - increase width of distal lobe promotion changed from 1.2 to 1.8  Eqn 49
                .* pro (0.1, id_cyc_l  .* pro(1.5,id_lip_l))... % increase width of dorsal petal    Eqn 44
                ...%%%%%% New added equations for Rebocho et al, for div mutant%%%%%
                .* pro (12, id_late_p .*s_rim_p .* pro (5, id_lpdp_l.* inh (100, id_lpdp_l > 0.644).* pro (5, id_hinge_l)).* inh (10,s_rad_l>0.028).*inh (100, (s_rad_l>0.5) .* inh (100, id_dich_l)))...% horizontal arms of orthogonal anisotropy
                .* inh (2, id_late_p .*id_lat_p .* (id_plt_p +id_lip_p).* pro (1, s_rim_l).*inh (100, id_rad_l).* inh (100, id_hinge_l))...%Keep petal junction regions narrow
                .* inh (2, id_late_p .* id_med_l.* (id_plt_l +id_lip_l).* pro (1, s_rim_l).*inh (100, id_rad_l).* inh (100, id_hinge_l))...%Keep MED regions narrow
                .* inh (4, s_rad_l .* id_hinge_l.* inh (100, id_tube_l .* inh (100, id_plt_l)).* pro (2, s_rim_l).* inh (10, id_prox_l).* inh (100, id_dtl_l .* inh (100, id_lip_l)))... % keeping the hinge narrow
                ...%%%%%% New added equations for Rebocho et al, for wild-type%%%%%
                .* pro (6, id_late_p .* id_div_l.* s_secvein_p .* id_lip_l.* inh (100, s_rim_l.^2))...%extend horizontal arms of orthogonal anisotropy to the LIP region
                .* inh (1, id_late_p.*(id_div_l>0.97) .* pro (2, s_rim_l>0.2).* inh (1, id_tube_l).* inh (10, id_lipcliff_l).* inh (100, id_dtl_l .* inh (100, id_lip_l)))... %keep ventral petal narrow but in the lipcliff regio
                .* inh (1, id_late_p .* id_mlobe_l .* (id_div_l < 0.995 & id_div_l >0.7));% help make lateral lip cheek protrude
            ...%%%%%%equations removed from published Fig 9B model7
                ....* inh (0.3, id_late_p .* id_lobe_l .* id_med_l)... %keep medial lobe narrow   Eqn 50
                ....* inh (1.3, id_late_p .* id_div_l .* inh(2.5, id_lobe_l.* inh (2,id_lpb_l)).* pro(1, id_plt_l))...  %narrow ventral petals    Eqn 51 %%%%%%% FINE TUNING %%%%%%
                ...%.* inh(6, id_early_p .* id_lts_l.*id_med_l )... %  keep laterals narrow initially   Eqn 45
                ...%.* pro (0.2, id_late_p .* id_lts_l .* id_med_l .*  id_plt_l); %promote growth of lateral width   Eqn 52   ADDD Fig 6G Model 3
                
            % compute the growth rates
            kapar=Kpar ...% Eqn 300
                .* inh (1, id_dich_l.* id_lat_l  .* id_lpb_l)...  % bend dorsal petals Eqn 33
                .* inh (1, id_early_p .*(id_cyc_l + id_div_l) .*id_mlobe_l .* inh (5, id_dich_l))...  %bend up distal lobes   Eqn 34
                .* inh (0.2, id_div_l .* id_lpb_l)... % Modified for Rebocho et al - bend out lips throughout all stages and uniformly MODIFIED FROM - .* inh (0.2, id_div_l .* id_lpb_l.* inh (5, id_med_l))  Eqn 35
                +0.04* id_early_p .* (s_rim_l>0.2) .* inh (2, id_cyc_l .* inh (25, id_dich_l)).* inh (1, id_hinge_l); % Modified for Rebocho et al - bend lower petals  - MODIFIED from + 0.05* id_early_p .* id_rim_l .* ( 0.8*id_div_l + 0* id_hinge_l.* id_med_l  + 0.5*(id_cyc_l +id_dich_l)); % Eqn 301
            kbpar=Kpar ... % Eqn 310
                .* pro (2, id_dich_l .* id_lat_l  .* id_lpb_l)... %Eqn 37 promote abaxial growth of lip triangle in between dorsal petal
                .* inh (0.5, (s_rim_l>0.15) .* pro (1, id_div_l).*inh (100, id_rad_l) .* inh (10, s_rad_l)); % NEW Rebocho et al - lower petal furrow to allow arching outwards
            kaper=Kper ... % Eqn 320
                .* pro (1, id_early_p .* id_div_l .* id_lpb_l)...%Modified for Rebocho et al -bend out lips - MIDIDIED FROM  Eqn 47 .* pro (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l))
                .* inh (0.5, id_early_p.*(id_div_l>0.97).*id_lipcliff_l);% slightly inflate the lower petal lip
            kbper=Kper ... % Eqn 330
                .* inh (1, id_early_p .* id_div_l .* id_lpb_l)...%Modified for Rebocho et al -bend out lips - MIDIDIED FROM  Eqn 48 .* inh (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l))
                .* inh (2, id_late_p.*(id_div_l>0.97).*id_lip_l); %slightly inflate the ventral petal lip
            knor =0.003;
            v_kaniso_p = log( kbpar_p./kbper_p);% visualise growth anisotropy in ln scale
            v_karea_p = kbper_p + kbpar_p;% visualise areal growth rate
            
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif strcmpi(modelnumber,'Figure 10K Wild type corolla') || strcmpi(modelnumber,'Fig 9B Model 7M2') || strcmpi(modelnumber,'Fig 9B Model 7M3')
        % P EQUATIONS
        % START POLARISER section (it is slightly mixed with gene networks)
        m = leaf_mgen_conductivity( m, 's_distorg', 0.002); %was 0.1 distorg_diff);% diffusion constant
        m = leaf_mgen_absorption( m, 's_distorg', 0.05);  % it will not decay everywhere
        m.globalProps.mingradient=0.01; % this ensures that all polariser gradients are frozen
        m = leaf_mgen_conductivity( m, 'Polariser', 0.002);% diffusion constant WAS 0.001
        m = leaf_mgen_dilution( m, 'Polariser', false );% it will not dilute with growth
        m = leaf_mgen_absorption( m, 'Polariser', 0.1);  % it will not decay everywhere
        
        m = leaf_setproperty (m, ...
            'twosidedpolarisation', true, ... % the mesh has two-sided polarisation
            'mingradient', 0.0, ... % No threshold for freezing the polariser
            'usedpolfreezebc', true); % a rather subtle choice of a feature of the behaviour of frozen gradients. It probably doesn't make much difference.
        
        m = leaf_setpolfrozen (m, false); %Initially the gradient is not frozen anywhere
        
        %\latex \subsection{PRN (changes)}
        % POLARISER PRODUCTION
        if realtime>230 % From time 230  generate polariser
            id_proxorg_p = (id_prox_l >0.8);                                             % Eqn 2  CHECKED
            id_proxorg_l = id_proxorg_p * id_proxorg_a;
        end
        
        % POLARISER DECAY
        % switch on distorg with id_dist in the absence of id_lat
        if abs(realtime - 230) <0.5*dt % At time 230 hrs
            id_distorg_p = 0.25* id_dist_l .* inh (100, id_lat_l)+ 0.1*id_dich_l .* id_dist_l .* id_cyc_l;           % Eqn 3 (used in production equation)
            id_distorg_l = id_distorg_p * id_distorg_a;
            m.userdata.all_regions_ready=true;
            v_dichdistorg_p=v_dichdistorg_p.*id_dich_l.* id_dist_l .* id_cyc_l;
        end
        
        if abs(realtime - 240) <0.5*dt %At 240 hrs
            
            %S_RAD NEW Rebocho et al%
            s_rad_p(:) = 0;
            s_rad_p(:) = 2*id_rad_p (:);
            m = leaf_mgen_conductivity( m, 'S_RAD', 0.002 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0.1 );
        end
        
        m.mgen_production(:,polariser_i) = 0.1 * (id_proxorg_l + 1 - P .* id_distorg_l); % Eqn 19 decay part of WAS S_DISTORG BY MISTAKE
        
        % LATER MODULATION OF DISTORG AND DIST
        % s_dist is activated by id_dist
        m.mgen_production(:,s_dist_i) = 0.1 * id_dist_l.* inh (100, id_lat_l);          % Eqn 17  CHECKED
        
        % s_distorg is activated by id_distorg
        m.mgen_production(:,s_distorg_i) = 0.1 * id_distorg_l ;                         % Eqn 15 (but  CHECKED
        
        if abs(realtime - 340) <0.5*dt % At time 340 hrs
            id_distorg_p = s_distorg_l;                                                 % Eqn 12  CHECKED
            id_distorg_l = id_distorg_p .* id_distorg_a;
            id_dist_p = s_dist_l;                                                       % Eqn 13  CHECKED (NOT NEEDED IN THIS MODEL)
            id_dist_l = id_dist_p .* id_dist_a;
        end
        % End Polariser
        %\latex \subsection{GRN (changes)}
        
        % id_rad is activated by id_cyc and id_dich at ALL TIMES
        id_rad_p = id_cyc_l +  0.3* id_dich_l;                                           % Eqn 1 % CHECKED
        id_rad_l = id_rad_p * id_rad_a;
        
        if abs(realtime - 250) <0.5*dt % At time 250 hrs
            %id_cyc and id_dich are inhibited by id_div in absence of id_rad.
            %This is to ensure that in a id_rad mutant,
            %id_cyc and id_dich are inactivated. (Also see constitutive cyc
            %below)
            id_cyc_p = id_cyc_l.* inh (4, id_div_l );                                    % Eqn 7 % CHECKED
            id_cyc_l = id_cyc_p * id_cyc_a;
            
            id_dich_p = id_dich_l .* inh (4, id_div_l );                                 % Eqn 8 % CHECKED
            id_dich_l = id_dich_p * id_dich_a;
            
            %S_RAD gradient stabilises- NEW Rebocho et al%
            m = leaf_mgen_conductivity( m, 'S_RAD', 0 );
            m = leaf_mgen_absorption( m, 'S_RAD', 0 );
            
            %switch on id_div_NEW Rebocho et al%
            id_div_p=zeros(size(id_div_p));
            id_div_p=inh(100, id_rad_l).* inh (10, s_rad_p);                           % Eqn 5 % CHECKED
            id_div_l=id_div_p.*id_div_a;
            
            %Specification of hinge region_NEW Rebocho et al%
            id_hinge_p = (s_rad_l<0.39).* (s_rad_l>0.03);
            id_hinge_l = id_hinge_p .* id_hinge_a;
            
            %Specification of LP_DP region_NEW Rebocho et al%
            id_lpdp_p = (s_rad_l<0.55).* (s_rad_l>0.2).* id_lat_l;
            id_lpdp_l = id_lpdp_p .* id_lpdp_a;
            
            %id_msinus_NEW Rebocho et al%
            id_msinus_p = (id_lat_l>0.65).* id_lobe_l.* inh (10, s_rim_l) ; %graded expression of ID MSINUS from sinus up to rim
            id_msinus_l = id_msinus_p .* id_msinus_a;
            
            id_lipvisual_p = id_lip_l.*inh(100, id_rim_l);
            
        end
        
        if realtime >= 250 && realtime<352.5 %Freeze polarity in the tube region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, id_tube_p, id_tube_p);
        end
        
        if realtime >= 340 % Reorientation of polarity at the ventral-lateral sinuses_NEW Rebocho et al%
            P(id_msinus_p.* id_div_l >0.7) = 0.9;
            m.morphogenclamp((id_msinus_p.* id_div_l> 0.7), polariser_i) = 1;
            m = leaf_setpolfrozen (m, false, id_lip_p);
        end
        
        if realtime >= 352.5 %Freeze polarity in the tube and Lip region_NEW Rebocho et al%
            m = leaf_setpolfrozen (m, (id_lip_p + id_tube_l), (id_lip_p + id_tube_l));
        end
        
        % End Gene Networks
        
        %\latex \subsection{KRN (changes)}
        if realtime >= 250 && realtime<570  % Between time 250 570
            % Compute growth rates
            
            % Note: the factors are vectors, one element per vertex (node)
            % as a result, with each statement, different things happen in
            % different regions of the model.
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % need to delete effect of cenorg on polarity
            
            Kpar=... % Eqn 260
                0.013 ... % k1.1 Background growth parallel to the polariser
                ...%%%%%equations maintained from snapdragon published Fig 9B model7%%%%%%
                .* inh (0.2, id_prox_l)... %k1.2 keeps base small
                .* pro (0.4,id_tube_l .* inh (100, id_plt_l))... %k1.3 promote growth of tube
                .* pro (1.4, id_early_p .* id_lip_l .* id_dich_l.* id_lat_l)... %k1.4b Modified for Rebocho et al - To grow the dorsal lip a bit more
                .* inh (0.5, id_late_p .* id_med_l.* id_lobe_l.* id_cyc_l)... %k1.6a Modified for Rebocho et al - to separate the control of the growth of dorsal lobe from the lower petal lobes
                .* inh (0.3, id_late_p .* id_lobe_l.*pro (0.5, id_med_l).* inh (5, id_dist_l ).* inh (5,id_lipcliff_l).* inh (100, id_cyc_l))... %k1.6b NEW Rebocho et al - TO GROW VENTRAL LOBE MORE
                .* inh (1, id_rad_l .*id_lip_l .* inh(6, id_dich_l .* id_lat_l))... %k1.8 Modified for Rebocho et al - prevent lip growing at edge of dorsals
                .* inh (1, id_rad_l .*id_plt_l .* inh (15, id_dich_l .* inh (5, id_lat_l) ) .* inh (30, id_cyc_l .* inh (40, id_lat_l.^2)))... %k1.9 prevent palate growing at edge of dorsals   Eqn 27 % CHECKED
                .* pro (0.3, id_cyc_l .* id_dtl_l  .* inh (0.5, id_dich_l))... %k1.10 promote growth of lateral dorsal lobe         Eqn 28 % CHECKED
                .* pro (0.45, (id_cyc_l + 0.2* id_dich_l) .* id_plt_l)... %k1.11 promote medial palate growth with CYC              Eqn 29 % CHECKED
                ...%%%%%% New added equations for Rebocho et al
                .* pro (2, id_early_p .* id_lipdistal_l.*inh (5, s_rad_l) .*inh (100, id_rad_l).* inh (0.7, id_div_l))...%k1.16a grow lipdistal only at ealy stages
                .* pro (0.5, id_late_p .* id_lipdistal_l .*inh (5, s_rad_l).*inh (100, id_rad_l).* inh (0.7, id_div_l))...%k1.17a grow lipdistal only at late stages
                .* inh (2, s_rad_l .*(id_lip_p + 0.5*id_uptube_l).* inh (100,s_rad_l > 0.25))...%k1.18 Keep hinge region small
                .* pro (0.35, id_late_p .* id_dtl_p .* id_hinge_l .*inh (100, id_lip_l))... %k1.19 grow distal lobes in hinge region
                .* pro (0.3, id_early_p .* id_lobe_l.* inh (100, id_lip_l).*inh (1, id_rad_l))...%k1.20 Grow lobe region at early stages - necessary to ovatined bigger lobes in all petals
                .* inh (0.5, id_late_p.* id_med_l .* (id_lip_l + id_uptube_l) .* pro (2, s_rim_l) .*inh (100, id_rad_l))... %k1.21 keep midvein region narrow
                .* pro (0.8, id_late_p.* id_lat_l .* (id_lip_l + id_uptube_l ) .*inh (100, s_rad_l>0.4).* pro(1, id_lpdp_l).* inh (100, id_lip_l.*id_lpdp_l))...%k1.22 vertical arms of orthogonal anisotropy
                .* inh (2, id_late_p .*(s_rim_l.^2).*inh (100, id_rad_l))... %k1.23 Keep rim region short
                .* pro (2.5, id_early_p.* id_div_l.*id_lipcliff_p .*inh (100, id_rim_l))...%k1.24 grow lipcliff only at early stages
                .* pro (0.2, id_late_p.* id_div_l.*id_lipcliff_p .*inh (100, id_rim_l))....%k1.25 grow lipcliff only at late stages ...% GROW PLT EARLY IS ESSENTIAL FOR BENDING THEPETALS AND HELPING WITH SHAPE
                .* pro (2, id_early_p .* id_div_l .* id_plt_l )... .* pro (0.4, id_uptube_p))... %k1.26 promote growth in palate with DIV
                .* pro (0.45, id_late_p .* id_div_l .* id_plt_l .* inh (2, id_uptube_l .* inh (100, id_lpdp_l)))... %k1.27 promote growth in palate with DIV                
                .* inh (0.8, id_late_p.* id_med_l .* (0.5*id_plt_p + id_lipcliff_l).*(id_div_l>0.97)); %k1.28 make ventral cleft

            Kper=... % Eqn 250
                0.0075 ... % Background growth normal to the polariser   Emqn 201
                ...%%%%%equations maintained from snapdragon published Fig 9B odel7%%%%%%
                .* inh (0.2, id_prox_l)... %keep base small   Eqn 43
                .* pro (2, id_early_p .* id_dist_l .* inh (20, id_lat_l))... %  widen distal lobe    Eqn 46
                .* pro (1, id_late_p .* id_dist_l .* pro (1.8, (id_cyc_l + id_div_l)))... % Modified for Rebocho et al - increase width of distal lobe promotion changed from 1.2 to 1.8  Eqn 49
                .* pro (0.1, id_cyc_l  .* pro(1.5,id_lip_l))... % increase width of dorsal petal    Eqn 44
                ...%%%%%% New added equations for Rebocho et al, for div mutant%%%%%
                .* pro (12, id_late_p .*s_rim_p .* pro (5, id_lpdp_l.* inh (100, id_lpdp_l > 0.644).* pro (5, id_hinge_l)).* inh (10,s_rad_l>0.028).*inh (100, (s_rad_l>0.5) .* inh (100, id_dich_l)))...% horizontal arms of orthogonal anisotropy
                .* inh (2, id_late_p .*id_lat_p .* (id_plt_p +id_lip_p).* pro (1, s_rim_l).*inh (100, id_rad_l).* inh (100, id_hinge_l))...%Keep petal junction regions narrow
                .* inh (2, id_late_p .* id_med_l.* (id_plt_l +id_lip_l).* pro (1, s_rim_l).*inh (100, id_rad_l).* inh (100, id_hinge_l))...%Keep MED regions narrow
                .* inh (4, s_rad_l .* id_hinge_l.* inh (100, id_tube_l .* inh (100, id_plt_l)).* pro (2, s_rim_l).* inh (10, id_prox_l).* inh (100, id_dtl_l .* inh (100, id_lip_l)))... % keeping the hinge narrow
                ...%%%%%% New added equations for Rebocho et al, for wild-type%%%%%
                .* pro (6, id_late_p .* id_div_l.* s_secvein_p .* id_lip_l.* inh (100, s_rim_l.^2))...%extend horizontal arms of orthogonal anisotropy to the LIP region
                .* inh (1, id_late_p.*(id_div_l>0.97) .* pro (2, s_rim_l>0.2).* inh (1, id_tube_l).* inh (10, id_lipcliff_l).* inh (100, id_dtl_l .* inh (100, id_lip_l)))... %keep ventral petal narrow but in the lipcliff regio
                .* inh (1, id_late_p .* id_mlobe_l .* (id_div_l < 0.995 & id_div_l >0.7));% help make lateral lip cheek protrude
            ...%%%%%%equations removed from published Fig 9B model7
                ....* inh (0.3, id_late_p .* id_lobe_l .* id_med_l)... %keep medial lobe narrow   Eqn 50
                ....* inh (1.3, id_late_p .* id_div_l .* inh(2.5, id_lobe_l.* inh (2,id_lpb_l)).* pro(1, id_plt_l))...  %narrow ventral petals    Eqn 51 %%%%%%% FINE TUNING %%%%%%
                ...%.* inh(6, id_early_p .* id_lts_l.*id_med_l )... %  keep laterals narrow initially   Eqn 45
                ...%.* pro (0.2, id_late_p .* id_lts_l .* id_med_l .*  id_plt_l); %promote growth of lateral width   Eqn 52   ADDD Fig 6G Model 3
                
            % compute the growth rates
            kapar=Kpar ...% Eqn 300
                .* inh (1, id_dich_l.* id_lat_l  .* id_lpb_l)...  % bend dorsal petals Eqn 33
                .* inh (1, id_early_p .*(id_cyc_l + id_div_l) .*id_mlobe_l .* inh (5, id_dich_l))...  %bend up distal lobes   Eqn 34
                .* inh (0.2, id_div_l .* id_lpb_l)... % Modified for Rebocho et al - bend out lips throughout all stages and uniformly MODIFIED FROM - .* inh (0.2, id_div_l .* id_lpb_l.* inh (5, id_med_l))  Eqn 35
                +0.04* id_early_p .* (s_rim_l>0.2) .* inh (2, id_cyc_l .* inh (25, id_dich_l)).* inh (1, id_hinge_l); % Modified for Rebocho et al - bend lower petals  - MODIFIED from + 0.05* id_early_p .* id_rim_l .* ( 0.8*id_div_l + 0* id_hinge_l.* id_med_l  + 0.5*(id_cyc_l +id_dich_l)); % Eqn 301
            kbpar=Kpar ... % Eqn 310
                .* pro (2, id_dich_l .* id_lat_l  .* id_lpb_l)... %Eqn 37 promote abaxial growth of lip triangle in between dorsal petal
                .* inh (0.5, (s_rim_l>0.15) .* pro (1, id_div_l).*inh (100, id_rad_l) .* inh (10, s_rad_l)); % NEW Rebocho et al - lower petal furrow to allow arching outwards
            kaper=Kper ... % Eqn 320
                .* pro (1, id_early_p .* id_div_l .* id_lpb_l)...%Modified for Rebocho et al -bend out lips - MIDIDIED FROM  Eqn 47 .* pro (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l))
                .* inh (0.5, id_early_p.*(id_div_l>0.97).*id_lipcliff_l);% slightly inflate the lower petal lip
            kbper=Kper ... % Eqn 330
                .* inh (1, id_early_p .* id_div_l .* id_lpb_l)...%Modified for Rebocho et al -bend out lips - MIDIDIED FROM  Eqn 48 .* inh (1, id_early_p .* id_div_l .* id_lpb_l.* inh (3, id_med_l))
                .* inh (2, id_late_p.*(id_div_l>0.97).*id_lip_l); %slightly inflate the ventral petal lip
            knor =0.003;
            v_kaniso_p = log( kbpar_p./kbper_p);% visualise growth anisotropy in ln scale
            v_karea_p = kbper_p + kbpar_p;% visualise areal growth rate
            
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if realtime >= 250 && realtime<570  % Between time 250 570
        kapar_p=kapar;
        kbpar_p=kbpar;
        kaper_p=kaper;
        kbper_p=kbper;
        knor_p =knor;
    end
    % Growth factors cannot be negative.
    if realtime>250
        kapar_p = max( kapar_p, 0 );
        kbpar_p = max( kbpar_p, 0 );
        kaper_p = max( kaper_p, 0 );
        kbper_p = max( kbper_p, 0 );
    end
    % Reconstruct the old growth and bend morphogens.
    kpar_p = (kapar_p + kbpar_p)/2;
    kper_p = (kaper_p + kbper_p)/2;
    bendpar_p = kbpar_p - kpar_p;
    bendper_p = kbper_p - kper_p;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    if realtime>570 % finished growth
        % just increment time to allow flattened results to be saved
        kapar_p = 0;
        kbpar_p = 0;
        kaper_p = 0;
        kbper_p = 0;
    end
    
    % dissect the corolla for flattening
    %         if realtime==570 %(570>realtime-dt) && (570<realtime+dt)
    %                     lateral_protect=double(id_div_p>0.2 & id_div_p<0.536);%536);
    %                     f_seam4_p=lateral_protect;
    %                     dorsal_protect=(id_div_p<0.01);
    %                     f_seam4_p=dorsal_protect;
    %                     f_seam_p=f_seppetals_p + f_sepbase_p + f_seplobes2_p .* (f_seam3_p+lateral_protect);
    %                     m=leaf_set_seams(m,f_seam_p);
    %         elseif realtime==(570+dt)
    % USE this for separate parts
    f_seam_p(:)=0;
    f_seam_p=f_seppetals_p + f_sepbase_p + f_seplobes2_p ;
    %f_seam_p=f_seppetals_p + f_sepbase_p ;
    m=leaf_set_seams(m,f_seam_p);
    %         elseif realtime==(570+2*dt)
    %             f_seam_p=f_seppetals_p + f_sepbase_p + f_seplobes2_p ;
    %             m=leaf_set_seams(m,f_seam_p);
    %         end
    
    % Calculate the area, excluding the base.
    areaCells = any( id_prox_p( m.tricellvxs )==0, 2 );
    % These are the finite elements which have at least one vertex where
    % the id_prox morphogen is zero.
    
    nonBaseArea = sum( m.cellareas( areaCells ) );
    % This is the total area of these finite elements.
    
    fprintf( 1, 'Snapdragon area At %.3f, excluding the base = %.3f sq mm, including base = %.3f sq mm\n', ...
        realtime, nonBaseArea, sum( m.cellareas ) );
    
    % Calculate the area, excluding the base.
    nonBaseCells = any( id_prox_p( m.tricellvxs )==0, 2 );
    % These are the finite elements which have at least one vertex where
    % the id_prox morphogen is zero.
    nonBaseArea = sum( m.cellareas( nonBaseCells ) );
    % This is the total area of these finite elements.
    m.userdata.areas(m.globalDynamicProps.currentIter+1) = nonBaseArea;
    m.userdata.times(m.globalDynamicProps.currentIter+1) = realtime;
    
end