cell2cell_transport/SmithPhyllotaxis.cpp

Go to the documentation of this file.

/*
 * Copyright 2011-2016 Universiteit Antwerpen
 *
 * Licensed under the EUPL, Version 1.1 or  as soon they will be approved by
 * the European Commission - subsequent versions of the EUPL (the "Licence");
 * You may not use this work except in compliance with the Licence.
 * You may obtain a copy of the Licence at: http://ec.europa.eu/idabc/eupl5
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the Licence is distributed on an "AS IS" basis,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the Licence for the specific language governing
 * permissions and limitations under the Licence.
 */
#include "SmithPhyllotaxis.h"

#include "bio/Cell.h"
#include "bio/Wall.h"
#include "bio/ReduceCellWalls.h"

namespace SimPT_Default {
namespace CellToCellTransport {

using namespace std;
using namespace boost::property_tree;

SmithPhyllotaxis::SmithPhyllotaxis(const CoreData& cd)
{
        Initialize(cd);
}

void SmithPhyllotaxis::Initialize(const CoreData& cd)
{
        // Keep a reference to core data
        m_cd = cd;

        // Assign parameters from core data to this model
        const auto& p = m_cd.m_parameters->get_child("smith_phyllotaxis");
        m_D = p.get<double>("D");
        m_T = p.get<double>("T");
        m_k_T = p.get<double>("k_T");
        m_c = p.get<double>("c");
}

void SmithPhyllotaxis::operator()(Wall* w, double* dchem_c1, double* dchem_c2)
{
        // Since we're looping over walls and not over cells special care needs to
        // be taken to ensure the cell indices are consistent with cell-based
        // equations "on paper".
        //
        // By convention, for the current wall, "cell 1" is i and "cell 2" is j.
        Cell * C_i = w->GetC1();
        Cell * C_j = w->GetC2();

        // There is no transport between a cell and the boundary polygon
        if (!(C_i->IsBoundaryPolygon() || C_j->IsBoundaryPolygon())) {
                // Shorthand notations for neighboring cell's chemicals ...
                const double IAA_i = C_i->GetChemical(0);
                const double IAA_j = C_j->GetChemical(0);
                const double PIN_i = C_i->GetChemical(1);
                const double PIN_j = C_j->GetChemical(1);
                // ... and for other common factors
                const double l_ij = w->GetLength();
                const double V_i = C_i->GetArea();
                const double V_j = C_j->GetArea();

                // IAA diffusion (passive flux) between neighboring cells
                // (Note that: diff_i_to_j = -diff_j_to_i)
                const double diff_j_to_i = l_ij * m_D * (IAA_j - IAA_i);
                dchem_c1[0] += diff_j_to_i / V_i; // dIAA/dt in cell i
                dchem_c2[0] -= diff_j_to_i / V_j; // dIAA/dt in cell j

                // Returns one term of: \sum_{b \in N(a)} l_ab * exp(c * IAA_b)
                // When used as func in ReduceCellWalls, the entire sum is calculated.
                function<double(Cell *, Cell *, Wall *)> P_ab_denom_term = \
                        [this](Cell * C_a, Cell * C_b, Wall * W_ab)->double {
                        // Make sure neither of the cells is the boundary polygon
                        // (technically speaking C_a won't be the boundary since it's
                        // already checked above. Double checking can do no harm though)
                        if (!(C_a->IsBoundaryPolygon() || C_b->IsBoundaryPolygon())) {
                                const double IAA_b = C_b->GetChemical(0);
                                const double l_ab = W_ab->GetLength();
                                return l_ab * exp(m_c * IAA_b);
                        } else {
                                return 0.0;
                        }
                };

                // IAA active transport (PIN1 mediated) over the current wall
                // (Note that: actr_i_to_j = -actr_j_to_i)
                const double P_ij = PIN_i * l_ij * exp(m_c * IAA_j) / ReduceCellWalls(C_i, P_ab_denom_term);
                const double P_ji = PIN_j * l_ij * exp(m_c * IAA_i) / ReduceCellWalls(C_j, P_ab_denom_term);
                const double actr_j_to_i = P_ji * IAA_j * IAA_j / (1.0 + m_k_T * IAA_i * IAA_i) \
                        - P_ij * IAA_i * IAA_i / (1.0 + m_k_T * IAA_j * IAA_j);
                dchem_c1[0] += m_T * actr_j_to_i / V_i; // dIAA/dt in cell i
                dchem_c2[0] -= m_T * actr_j_to_i / V_j; // dIAA/dt in cell j
        }
}

} // namespace
} // namespace