cell_housekeep/Meinhardt.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 "Meinhardt.h"

#include "bio/Cell.h"
#include "bio/Mesh.h"

#include <cmath>

namespace SimPT_Default {
namespace CellHousekeep {

using namespace std;
using namespace boost::property_tree;

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

void Meinhardt::Initialize(const CoreData& cd)
{
        m_cd = cd;
        auto& p = m_cd.m_parameters;

        m_cell_base_area               = p->get<double>("cell_mechanics.base_area");
        m_cell_expansion_rate          = p->get<double>("cell_mechanics.cell_expansion_rate");
        m_constituous_expansion_limit  = p->get<unsigned int>("meinhardt.constituous_expansion_limit");
        m_division_ratio               = p->get<double>("cell_mechanics.division_ratio");
        m_elastic_modulus              = p->get<double>("cell_mechanics.elastic_modulus");
        m_response_time                = p->get<double>("cell_mechanics.response_time");
        m_time_step                    = p->get<double>("model.time_step");
        m_vessel_expansion_rate        = p->get<double>("meinhardt.vessel_expansion_rate");
        m_vessel_inh_level             = p->get<double>("meinhardt.vessel_inh_level");
        m_viscosity_const              = p->get<double>("cell_mechanics.viscosity_const");

        m_area_incr                    = m_cell_expansion_rate * m_time_step;
        m_div_area                     = m_division_ratio * m_cell_base_area;
}

void Meinhardt::operator()(Cell* cell)
{
        const string ham_select = m_cd.m_parameters->get<string>("model.mc_hamiltonian");

        const double chem_0            = cell->GetChemical(0);
        const double chem_3            = cell->GetChemical(3);
        const double level             = m_vessel_inh_level;
        const unsigned int limit       = m_constituous_expansion_limit;
        const double t_area            = cell->GetTargetArea();
        const double t_length          = cell->GetTargetLength();
        const unsigned int cell_count  = m_cd.m_mesh->GetCells().size();

        // Cell expansion is inhibited by substrate (chem 3).
        double update_t_area = 0.0;
        if (!limit || cell_count < limit) {
                update_t_area = t_area + m_area_incr;
        } else {
                if (chem_0 < 0.5) {
                        const double incr = max((1.0 - level * chem_3) * m_area_incr, 0.0);
                        update_t_area = t_area + incr;
                } else {
                        update_t_area = t_area + m_vessel_expansion_rate;
                }
        }
        const double update_t_length  = t_length * sqrt(update_t_area / t_area);

        cell->SetTargetArea(update_t_area);
        cell->SetTargetLength(update_t_length);

        if (ham_select == "ElasticWall") {
                // Update the rest length of wall
                // TODO: 1.50 (maximal extension) and 1.25 (irreversible extension) are hidden parameters
                for (list<Wall*>::iterator i = cell->GetWalls().begin(); i != cell->GetWalls().end(); ++i) {
                        Wall* w = *i;
                        if (w->GetLength() > 1.50 * w->GetRestLength()) {
                                // Irreversible extension
                                w->SetRestLength(w->GetLength() / 1.25);
                        }
                }
        }

        if (ham_select == "Maxwell") {
                // Update the solute of cell and the rest length of wall
                const double solute = m_viscosity_const * sqrt(cell->GetArea());
                const double update_solute = solute - (solute - cell->GetSolute()) * exp(-0.1 * m_time_step/m_response_time);
                cell->SetSolute(update_solute);

                for (list<Wall*>::iterator i = cell->GetWalls().begin(); i != cell->GetWalls().end(); ++i) {
                        Wall* w = *i;

                        const double update_rest_length = w->GetLength()
                                - (w->GetLength() - w->GetRestLength())
                                        * exp(-0.01 * m_elastic_modulus * m_time_step / m_viscosity_const);
                        w->SetRestLength(update_rest_length);
                }
        }
}

} // namespace
} // namespace