diff --git a/src/build_models/leafplugin.cpp b/src/build_models/leafplugin.cpp
new file mode 100644
--- /dev/null
+++ b/src/build_models/leafplugin.cpp
@@ -0,0 +1,425 @@
+/*
+ *
+ * This file is part of the Virtual Leaf.
+ *
+ * The Virtual Leaf is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * The Virtual Leaf is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with the Virtual Leaf. If not, see .
+ *
+ * Copyright 2010 Roeland Merks.
+ *
+ */
+
+#include
+#include
+
+#include "simplugin.h"
+
+#include "parameter.h"
+
+#include "wallbase.h"
+#include "cellbase.h"
+#include "leafplugin.h"
+
+#include "far_mem_5.h"
+
+static const std::string _module_id("$Id$");
+
+bool batch = false;
+
+// To be executed after cell division
+void LeafPlugin::OnDivide(ParentInfo &parent_info, CellBase &daughter1, CellBase &daughter2) {
+ // PIN1 distributes between parent and daughter according to area
+ double area = daughter1.Area(), daughter_area = daughter2.Area();
+ double tot_area = area + daughter_area;
+
+ //chem[1]*=(area/tot_area);
+ //daughter.chem[1]*=(daughter_area/tot_area);
+
+ // For lack of detailed data, or a better rule, we assume that cells remain polarized
+ // after division
+
+ // So the PIN1 is redistributed according to the original polarization over the walls
+
+ // parent_info contains info about the parent
+ // redistribute the PIN in the endosome according to area
+
+ // "Fudge" rule: if one of the cells is at the boundary, remove all AUX1 in the other cell
+ if (daughter1.AtBoundaryP() && !daughter2.AtBoundaryP()) {
+ //daughter2.new_chem[2]=daughter2.chem[2]=0.;
+ daughter2.SetNewChem(2,0);
+ daughter2.SetChemical(2,0);
+ //daughter.new_chem[0]=daughter.chem[0]=0.;
+ //cerr << "Clearing daughter\n";
+ //for (list::const_iterator w=daughter.walls.begin();
+ // w!=daughter.walls.end();
+ // w++) {
+
+ // (*w)->setTransporter(&daughter, 1, 0.);
+
+ //}
+ //new_chem[2]=chem[2]=parent_info.PINendosome;
+ daughter1.SetNewChem(2,parent_info.PINendosome);
+ daughter1.SetChemical(2,parent_info.PINendosome);
+
+ } else {
+ if (daughter2.AtBoundaryP() && !daughter1.AtBoundaryP()) {
+
+ //new_chem[2]=chem[2]=0.;
+ daughter1.SetNewChem(2,0);
+ daughter1.SetChemical(2,0);
+
+ /*new_chem[0]=chem[0]=0.;
+ for (list::const_iterator w=walls.begin();
+ w!=walls.end();
+ w++) {
+
+ (*w)->setTransporter(this, 1, 0.);
+ }*/
+ //daughter2.chem[2]=parent_info.PINendosome;
+ daughter2.SetChemical(2,parent_info.PINendosome);
+ //cerr << "Clearing parent\n";
+
+ } else {
+ //daughter1.new_chem[2]=daughter1.chem[2] = parent_info.PINendosome*(area/tot_area);
+ daughter1.SetNewChem(2,parent_info.PINendosome*(area/tot_area));
+ daughter1.SetChemical(2, parent_info.PINendosome*(area/tot_area));
+ //daughter2.new_chem[2]=daughter2.chem[2] = parent_info.PINendosome*(daughter_area/tot_area);
+ daughter2.SetNewChem(2,parent_info.PINendosome*(daughter_area/tot_area));
+ daughter2.SetChemical(2,parent_info.PINendosome*(daughter_area/tot_area));
+
+ }
+ }
+
+ /*
+ // NB: Code commented out; not yet adapted to plugin format... RM 18/12/2009
+ // Now redistribute the membrane PINs according to the original polarization in the parent
+ // mmm... I'd like to have a better, biologically motivated rule for this,
+ // but for lack of something better... I hope I'm excused :-). Let's say the overall
+ // organization of the actin fibres is not completely destroyed after division...
+
+ // distribute wallPINs according to the circumference of the parent and daughter
+ double circ = Circumference( );
+ double daughter_circ = daughter.Circumference();
+ double tot_circ = circ + daughter_circ;
+
+ double wallPINs = (circ / tot_circ) * parent_info.PINmembrane;
+ double daughter_wallPINs = (daughter_circ / tot_circ) * parent_info.PINmembrane;
+
+ //cerr << "wallPINs = " << wallPINs << ", daughter_wallPINs = " << daughter_wallPINs << "sum = " << wallPINs + daughter_wallPINs << ", PINmembrane = " << parent_info.PINmembrane << endl;
+ // distrubute it according to the overall polarity
+ Vector polarization = parent_info.polarization.Normalised().Perp2D();
+
+ double sum=0.;
+ for (list::const_iterator w=walls.begin();
+ w!=walls.end();
+ w++) {
+
+ // distribute according to angle (0 degrees: maximum, 180 degrees minimum)
+ double tmp=InnerProduct((*w)->getWallVector(this),polarization); // move domain from [-1,1] to [0,1]
+
+ cerr << "[" << tmp << "]";
+ sum+=tmp;
+ //(*w)->setTransporter(this, 1,
+ }
+
+ //cerr << "Sum is " << sum << endl;
+ //double sum_wall_Pi = SumTransporters(1);
+
+ // After division, cells produce PIN1 (in intracellular storage) until total amount becomes Pi_tot
+ //SetChemical(1, par.Pi_tot - sum_wall_Pi );
+ //SetNewChem(1, Chemical(1));
+
+ //cerr << "[ " << sum_wall_Pi + Chemical(1) << "]";
+ */
+}
+
+void LeafPlugin::SetCellColor(CellBase &c, QColor &color) {
+
+ // Red: AUX1
+ // Green: Auxin
+ // Blue: van-3
+ // color.setRgb(chem[2]/(1+chem[2]) * 255.,(chem[0]/(1+chem[0]) * 255.),(chem[3]/(1+chem[3]) *255.) );
+ color.setRgb(c.Chemical(2)/(1+c.Chemical(2)) * 255.,(c.Chemical(0)/(1+c.Chemical(0)) * 255.),(c.Chemical(3)/(1+c.Chemical(3)) *255.) );
+
+
+}
+
+
+
+void LeafPlugin::CellHouseKeeping(CellBase &c) {
+
+ if (c.Boundary()==CellBase::None) {
+ if (c.Area() > par->rel_cell_div_threshold * c.BaseArea() ) {
+ //c.SetChemical(0,0);
+ c.Divide();
+ }
+
+ // expand if this is not a provascular cell
+ if (c.Chemical(3) < 0.7 ) {
+ c.EnlargeTargetArea(par->cell_expansion_rate);
+ }
+ }
+
+}
+
+void LeafPlugin::CelltoCellTransport(Wall *w, double *dchem_c1, double *dchem_c2) {
+
+ // leaf edge is const source of auxin
+ // (Neumann boundary condition: we specify the influx)
+ if (w->C2()->BoundaryPolP()) {
+ if (w->AuxinSource()) {
+ double aux_flux = par->leaf_tip_source * w->Length();
+ dchem_c1[0]+= aux_flux;
+
+ // dchem_c2 is undefined..!
+ return;
+ } else {
+ if (w->AuxinSink()) {
+
+ // efflux into Shoot Apical meristem
+ // we assume all PINs are directed towards shoot apical meristem
+ dchem_c1[0] -= par->sam_efflux * w->C1()->Chemical(0) / (par->ka + w->C1()->Chemical(0));
+
+ return;
+ } else {
+
+ // Active fluxes (PIN1 and AUX1 mediated transport)
+
+ // (Transporters measured in moles, here)
+ // efflux from cell 1 to cell 2
+ double trans12 = ( par->transport * w->Transporters1(1) * w->C1()->Chemical(0) / (par->ka + w->C1()->Chemical(0))
+ + par->aux1transport * w->C2()->Chemical(2) * w->C1()->Chemical(0) / (par->kaux1 + w->C1()->Chemical(0)) );
+
+ // efflux from cell 2 to cell 1
+ double trans21 = ( par->transport * w->Transporters2(1) * w->C2()->Chemical(0) / (par->ka + w->C2()->Chemical(0))
+ + par->aux1transport * w->C1()->Chemical(2) * w->C2()->Chemical(0) / (par->kaux1 + w->C2()->Chemical(0)) );
+
+ dchem_c1[0] += trans21 - trans12;
+ dchem_c2[0] += trans12 - trans21;
+ return;
+ }
+
+ }
+ }
+
+
+ if (w->C1()->BoundaryPolP()) {
+
+ if (w->AuxinSource()) {
+ double aux_flux = par->leaf_tip_source * w->Length();
+ dchem_c2[0] += aux_flux;
+ // dchem_c1 is undefined...!
+ return;
+ } else {
+
+ if (w->AuxinSink()) {
+
+
+ // efflux into Shoot Apical meristem
+ // we assume all PINs are directed towards shoot apical meristem
+
+ // no passive fluxes: outside is impermeable
+
+ // Active fluxes (PIN1 and AUX1 mediated transport)
+
+ // (Transporters measured in moles, here)
+ // efflux from cell 1 to cell 2
+ // assumption: no AUX1 in shoot apical meristem
+ double trans12 = ( par->transport * w->Transporters1(1) * w->C1()->Chemical(0) / (par->ka + w->C1()->Chemical(0)));
+ dchem_c1[0] += - trans12;
+
+ return;
+
+ //dchem_c2[0] -= par->sam_efflux * w->C2()->Chemical(0) / (par->ka + w->C2()->Chemical(0));
+
+ // return;
+ } else {
+
+ }
+ }
+ }
+
+
+ // Passive fluxes (Fick's law)
+ // only auxin flux now
+ // flux depends on edge length and concentration difference
+ for (int c=0;cLength() * ( par->D[c] ) * ( w->C2()->Chemical(c) - w->C1()->Chemical(c) );
+ dchem_c1[c] += phi;
+ dchem_c2[c] -= phi;
+ }
+ // Active fluxes (PIN1 and AUX1 mediated transport)
+
+ // (Transporters measured in moles, here)
+ // efflux from cell 1 to cell 2
+ double trans12 = ( par->transport * w->Transporters1(1) * w->C1()->Chemical(0) / (par->ka + w->C1()->Chemical(0))
+ + par->aux1transport * w->C2()->Chemical(2) * w->C1()->Chemical(0) / (par->kaux1 + w->C1()->Chemical(0)) );
+
+ // efflux from cell 2 to cell 1
+ double trans21 = ( par->transport * w->Transporters2(1) * w->C2()->Chemical(0) / (par->ka + w->C2()->Chemical(0))
+ + par->aux1transport * w->C1()->Chemical(2) * w->C2()->Chemical(0) / (par->kaux1 + w->C2()->Chemical(0)) );
+
+ dchem_c1[0] += trans21 - trans12;
+ dchem_c2[0] += trans12 - trans21;
+
+
+
+}
+void LeafPlugin::WallDynamics(Wall *w, double *dw1, double *dw2) {
+
+ // Cells polarize available PIN1 to Shoot Apical Meristem
+ if (w->C2()->BoundaryPolP()) {
+ if (w->AuxinSink()) {
+
+ dw1[0] = 0.; dw2[0] = 0.;
+ dw1[2] = 0.; dw2[2] = 0.;
+
+ // assume high auxin concentration in SAM, to convince PIN1 to polarize to it
+ // exocytosis regulated0
+ double nb_auxin = par->sam_auxin;
+ double receptor_level = nb_auxin * par->r / (par->kr + nb_auxin);
+
+ dw1[1] = par->k1 * w->C1()->Chemical(1) * receptor_level /( par->km + w->C1()->Chemical(1) ) - par->k2 * w->Transporters1(1);
+
+ dw2[1] = 0.;
+ return;
+
+ } else {
+ dw1[0]=dw2[0]=dw1[1]=dw2[1]=dw1[2]=dw2[2];
+ return;
+ }
+ }
+
+ if (w->C1()->BoundaryPolP()) {
+ if (w->AuxinSink()) {
+
+ dw1[0] = 0.; dw2[0] = 0.;
+ dw1[2] = 0.; dw2[2] = 0.;
+
+ // assume high auxin concentration in SAM, to convince PIN1 to polarize to it
+ // exocytosis regulated
+ double nb_auxin = par->sam_auxin;
+ double receptor_level = nb_auxin * par->r / (par->kr + nb_auxin);
+ dw2[1] = par->k1 * w->C2()->Chemical(1) * receptor_level /( par->km + w->C2()->Chemical(1) ) - par->k2 * w->Transporters2(1);
+
+ dw1[1] = 0.;
+ return;
+
+ } else {
+ dw1[0]=dw2[0]=dw1[1]=dw2[1]=dw1[2]=dw2[2];
+ return;
+ }
+ }
+
+
+
+ // PIN1 localization at wall 1
+ // Note: chemical 0 is Auxin (intracellular storage only)
+ // Chemical 1 is PIN1 (walls and intracellular storage)
+ //! \f$ \frac{d Pij/dt}{dt} = k_1 A_j \frac{P_i}{L_ij} - k_2 P_{ij} \f$
+ // Note that Pij is measured in term of concentration (mol/L)
+ // Pi in terms of quantity (mol)
+
+ double dPijdt1=0., dPijdt2=0.;
+
+ // normal cell
+ double auxin2 = w->C2()->Chemical(0);
+ double receptor_level1 = auxin2 * par->r / (par->kr + auxin2);
+
+ dPijdt1 =
+ // exocytosis regulated
+ par->k1 * w->C1()->Chemical(1) * receptor_level1 / ( par->km + w->C1()->Chemical(1) ) - par->k2 * w->Transporters1(1);
+
+ double auxin1 = w->C1()->Chemical(0);
+ double receptor_level2 = auxin1 * par->r / (par->kr + auxin1);
+
+ // normal cell
+ dPijdt2 =
+
+ // exocytosis regulated
+ par->k1 * w->C2()->Chemical(1) * receptor_level2 / ( par->km + w->C2()->Chemical(1) ) - par->k2 * w->Transporters2(1);
+
+ /* PIN1 of neighboring vascular cell inhibits PIN1 endocytosis */
+
+ dw1[0] = 0.; dw2[0] = 0.;
+ dw1[2] = 0.; dw2[2] = 0.;
+
+ dw1[1] = dPijdt1;
+ dw2[1] = dPijdt2;
+
+}
+
+double LeafPlugin::complex_PijAj(CellBase &here, CellBase &nb, Wall &w) {
+
+ // gives the amount of complex "auxinreceptor-Pin1" at the wall (at QSS)
+ //return here.Chemical(1) * nb.Chemical(0) / ( par->km + here.Chemical(1));
+
+ double nb_aux = (nb.BoundaryPolP() && w.AuxinSink()) ? par->sam_auxin : nb.Chemical(0);
+ double receptor_level = nb_aux * par->r / (par->kr + nb_aux);
+
+ return here.Chemical(1) * receptor_level / ( par->km + here.Chemical(1));
+
+}
+
+
+void LeafPlugin::CellDynamics(CellBase *c, double *dchem) {
+
+ double dPidt = 0.;
+
+ double sum_Pij = c->SumTransporters( 1 );
+
+ // exocytosis regulated:
+ // van3 expression reduces rate of PIN1 endocytosis
+ dPidt = -par->k1 * c->ReduceCellAndWalls( far_3_arg_mem_fun( *this, &LeafPlugin::complex_PijAj ) ) +
+ (c->Chemical(3) < 0.5 ? par->k2 : par->k2van3) * sum_Pij;
+
+ // production of PIN depends on auxin concentration
+ dPidt += (c->AtBoundaryP()?par->pin_prod_in_epidermis:par->pin_prod) * c->Chemical(0) - c->Chemical(1) * par->pin_breakdown;
+
+ /*if (c->AtBoundaryP()) {
+ dchem[2] = 0.01;
+ //cerr << "Making cell blue.\n";
+ } else {
+ dchem[2] = -0.1 * c->Chemical(2);
+ }*/
+
+ // no PIN production in SAM
+ if (c->Boundary() == CellBase::SAM) {
+ dchem[1]=0.;
+ dchem[0]= - par->sam_auxin_breakdown * c->Chemical(0);
+ dchem[2]=0.;
+ } else {
+
+ dchem[1] = dPidt;
+
+
+ // source of auxin
+ dchem[0] = par->aux_cons;
+
+ // auxin-induced AUX1 production, in the epidermis
+ dchem[2] = ( c->AtBoundaryP() ? par->aux1prod : par->aux1prodmeso ) * ( c->Chemical(0) / ( 1. + par->kap * c->Chemical(0) ) ) - ( par->aux1decay ) * c->Chemical(2) ;//: 0.;
+
+ // auxin-induced production of VAN-3? Autokatalysis?
+ //dchem[3] = par->van3prod * (c->Chemical(0) / (1. + par->kvp * c-> Chemical(0) ) )
+ double A = c->Chemical(0);
+ double van3 = c->Chemical(3);
+ dchem[3] = par->van3prod * A - par->van3autokat * van3 + van3*van3/(1 + par->van3sat * van3*van3 );
+ }
+}
+
+
+
+
+
+
+Q_EXPORT_PLUGIN2(leafplugin, LeafPlugin)