Skip to content
Snippets Groups Projects
Unverified Commit 397a66b5 authored by Sebastian Schmitt's avatar Sebastian Schmitt Committed by GitHub
Browse files

lmorpho (#1746)

Add the lmorpho utility back to Arbor

The original utility was removed because it was a maintenance burden with no users.
Now users are requesting the feature, this PR adds the old code, and updates it
to work with the new morphology description API.
parent c88b7be5
No related branches found
No related tags found
No related merge requests found
......@@ -556,3 +556,4 @@ install(
cmake/FindUnwind.cmake
DESTINATION "${cmake_config_dir}")
add_subdirectory(lmorpho)
add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp)
target_link_libraries(lmorpho PRIVATE arbor arbor-sup ext-tinyopt arborio)
# TODO: resolve public headers
target_link_libraries(lmorpho PRIVATE arbor-private-headers)
install(TARGETS lmorpho RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
# `lmporho` — generate random neuron morphologies
`lmorpho` implements an L-system geometry generator for the purpose of
generating a series of artificial neuron morphologies as described by a set of
stasticial parameters.
## Method
The basic algorithm used is that of [Burke (1992)][burke1992], with extensions
for embedding in 3-d taken from [Ascoli (2001)][ascoli2001].
A dendritic tree is represented by a piece-wise linear embedding of a rooted
tree into R³×R⁺, describing the position in space and non-negative radius.
Morphologies are represented as a sequence of these trees together with a point
and radius describing a (spherical) soma.
The generation process for a tree successively grows the terminal points by
some length ΔL each step, changing its orientation each step according to a
random distribution. After growing, there is a radius-dependent probability
that it will bifurcate or terminate. Parameters describing these various
probability distributions constitute the L-system model.
The distributions used in [Ascoli (2001)][ascoli2001] can be constant (i.e. not
random), uniform over some interval, truncated normal, or a mixture
distribution of these. In the present version of the code, mixture
distributions are not supported.
## Output
Generated morphologies can be output either in
[SWC format](http://research.mssm.edu/cnic/swc.html), or as 'parent vectors',
which describe the topology of the associated tree.
Entry _i_ of the (zero-indexed) parent vector gives the index of the proximal
unbranched section to which it connects. Somata have a parent index of -1.
Morphology information can be written one each to separate files or
concatenated. If concatenated, the parent vector indices will be shifted as
required to maintain consistency.
## Models
Two L-system parameter sets are provided as 'built-in'; if neither model is
selected, a simple, non-biological default model is used.
As mixture distributions are not yet supported in `lmorpho`, these models are
simplified approximations to those in the literature as described below.
### Motor neuron model
The 'motoneuron' model corresponds to the adult cat spinal alpha-motoneuron
model described in [Ascoli (2001)][ascoli2001], based in turn on the model of
[Burke (1992)][burke1992].
The model parameters in Ascoli are not complete; missing parameters were taken
from the corresponding ArborVitae model parameters in the same paper, and
somatic diameters were taken from [Cullheim (1987)][cullheim1987].
### Purkinke neuron model
The 'purkinke' model corresponds to the guinea pig Purkinje cell model from
[Ascoli (2001)][ascoli2001], which was fit to data from [Rapp
(1994)][rapp1994]. Somatic diameters are a simple fit to the three measurements
in the original source.
Some required parameters are not recorded in [Ascoli (2001)][ascoli2001];
the correlation component `diam_child_a` and the `length_step` ΔL are
taken from the corresponding L-Neuron parameter description in the
[L-Neuron database](l-neuron-database). These should be verified by
fitting against the digitized morphologies as found in the database.
Produced neurons from this model do not look especially realistic.
### Other models
There is not yet support for loading in arbitrary parameter sets, or for
translating or approximating the parameters that can be found in the
[L-Neuron database][l-neuron-database]. Support for parameter estimation
from a population of discretized morphologies would be useful.
## References
1. Ascoli, G. A., Krichmar, J. L., Scorcioni, R., Nasuto, S. J., Senft, S. L.
and Krichmar, G. L. (2001). Computer generation and quantitative morphometric
analysis of virtual neurons. _Anatomy and Embryology_ _204_(4), 283–301.
[DOI: 10.1007/s004290100201][ascoli2001]
2. Burke, R. E., Marks, W. B. and Ulfhake, B. (1992). A parsimonious description
of motoneuron dendritic morphology using computer simulation.
_The Journal of Neuroscience_ _12_(6), 2403–2416. [online][burke1992]
3. Cullheim, S., Fleshman, J. W., Glenn, L. L., and Burke, R. E. (1987).
Membrane area and dendritic structure in type-identified triceps surae alpha
motoneurons. _The Journal of Comparitive Neurology_ _255_(1), 68–81.
[DOI: 10.1002/cne.902550106][cullheim1987]
4. Rapp, M., Segev, I. and Yaom, Y. (1994). Physiology, morphology and detailed
passive models of guinea-pig cerebellar Purkinje cells.
_The Journal of Physiology_, _474_, 101–118.
[DOI: 10.1113/jphysiol.1994.sp020006][rapp1994].
[ascoli2001]: http://dx.doi.org/10.1007/s004290100201
"Ascoli et al. (2001). Computer generation […] of virtual neurons."
[burke1992]: http://www.jneurosci.org/content/12/6/2403
"Burke et al. (1992). A parsimonious description of motoneuron dendritic morphology […]"
[cullheim1987]: http://dx.doi.org/10.1002/cne.902550106
"Cullheim et al. (1987). Membrane area and dendritic structure in […] alpha motoneurons."
[rapp1994]: http://dx.doi.org/10.1113/jphysiol.1994.sp020006
"Rapp et al. (1994). Physiology, morphology […] of guinea-pig cerebellar Purkinje cells."
[l-neuron-database]: http://krasnow1.gmu.edu/cn3/L-Neuron/database/index.html
#include <algorithm>
#include <iostream>
#include <fstream>
#include <map>
#include <sstream>
#include <vector>
#include <tinyopt/tinyopt.h>
#include "lsystem.hpp"
#include "lsys_models.hpp"
#include <arbor/cable_cell.hpp>
#include <arborio/cableio.hpp>
#include <arborio/label_parse.hpp>
const char* usage_str =
"[OPTION]...\n"
"\n"
" -n, --count=N Number of morphologies to generate.\n"
" -m, --model=MODEL Use L-system MODEL for generation (see below).\n"
" --acc=FILE Output morphologies as ACC to FILE.\n"
" -h, --help Emit this message and exit.\n"
"\n"
"Generate artificial neuron morphologies based on L-system descriptions.\n"
"\n"
"If a FILE argument contains a '%', then one file will be written for\n"
"each generated morphology, with the '%' replaced by the index of the\n"
"morphology, starting from zero.\n"
"A FILE argument of '-' corresponds to standard output.\n"
"\n"
"Currently supported MODELs:\n"
" motoneuron Adult cat spinal alpha-motoneurons, based on models\n"
" and data in Burke 1992 and Ascoli 2001.\n"
" purkinje Guinea pig Purkinje cells, basd on models and data\n"
" from Rapp 1994 and Ascoli 2001.\n";
int main(int argc, char** argv) {
// options
int n_morph = 1;
std::optional<unsigned> rng_seed = 0;
lsys_param P;
bool help = false;
std::optional<std::string> acc_file;
std::pair<const char*, const lsys_param*> models[] = {
{"motoneuron", &alpha_motoneuron_lsys},
{"purkinje", &purkinje_lsys}
};
try {
for (auto arg = argv+1; *arg; ) {
bool ok = false;
ok |= (n_morph << to::parse<int>(arg, "--n", "--count")).has_value();
ok |= (rng_seed << to::parse<int>(arg, "--s", "--seed")).has_value();
ok |= (acc_file << to::parse<std::string>(arg, "-f", "--acc")).has_value();
const lsys_param* o;
if((o << to::parse<const lsys_param*>(arg, to::keywords(models), "-m", "--model")).has_value()) {
P = *o;
}
ok |= (help << to::parse(arg, "-h", "--help")).has_value();
if (!ok) throw to::option_error("unrecognized argument", *arg);
}
if (help) {
to::usage(argv[0], usage_str);
return 0;
}
std::minstd_rand g;
if (rng_seed) {
g.seed(rng_seed.value());
}
using namespace arborio::literals;
auto apical = apical_lsys;
auto basal = basal_lsys;
for (int i = 0; i < n_morph; ++i) {
const arb::segment_tree morph = generate_morphology(P.diam_soma, std::vector<lsys_param>{apical, basal}, g);
arb::label_dict labels;
labels.set("soma", "(tag 1)"_reg);
labels.set("basal", "(tag 3)"_reg);
labels.set("apical", "(tag 4)"_reg);
arb::decor decor;
arb::cable_cell cell(morph, labels, decor);
if(acc_file) {
std::string filename = acc_file.value();
const std::string toReplace("%");
auto pos = filename.find(toReplace);
if(pos != std::string::npos) {
filename.replace(pos, toReplace.length(), std::to_string(i));
}
if(filename == "-") {
arborio::write_component(std::cout, cell);
} else {
std::ofstream of(filename);
arborio::write_component(of, cell);
}
}
}
} catch (to::option_error& e) {
std::cerr << argv[0] << ": " << e.what() << "\n";
std::cerr << "Try '" << argv[0] << " --help' for more information.\n";
std::exit(2);
}
catch (std::exception& e) {
std::cerr << "caught exception: " << e.what() << "\n";
std::exit(1);
}
}
#include <cmath>
#include "lsystem.hpp"
#include "lsys_models.hpp"
static constexpr double inf = INFINITY;
// Predefined parameters for two classes of neurons. Numbers taken primarily
// from Ascoli et al. 2001, but some details (soma diameters for example)
// taken from earlier literature.
//
// Without mixture distribution support, these aren't entirely faithful to
// the original sources.
//
// Refer to `README.md` for references and more information.
// https://krasnow1.gmu.edu/cn3/L-Neuron/database/moto/bu/motoburk.prm
lsys_param make_alpha_motoneuron_lsys() {
lsys_param L;
// Soma diameter [µm].
L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
// Number of dendritic trees (rounded to nearest int).
L.n_tree = { 8.0, 16.0 }; // numtree
// Initial dendrite diameter [µm]. (Dstem)
L.diam_initial = { 3.0, 12.0 }; // initdiam
// Dendrite step length [µm]. (ΔL)
L.length_step = { 25 };
// Initial roll (intrinsic rotation about x-axis) [degrees].
L.roll_initial = { -180.0, 180.0 };
// Initial pitch (intrinsic rotation about y-axis) [degrees].
L.pitch_initial = { 0.0, 180.0 }; // biforient
// Tortuousness: roll within section over ΔL [degrees].
L.roll_section = { -180.0, 180.0 };
// Tortuousness: pitch within section over ΔL [degrees].
L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
// Taper rate: diameter decrease per unit length.
L.taper = { -1.25e-3 };
// Branching torsion: roll at branch point [degrees].
L.roll_at_branch = { 0.0, 180.0 };
// Branch angle between siblings [degrees].
L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
// Correlated child branch radius distribution parameters.
L.diam_child_a = -0.2087;
L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
// Termination probability parameters [1/µm].
L.pterm_k1 = 2.62e-2;
L.pterm_k2 = -2.955;
// Bifurcation probability parameters [1/μm].
L.pbranch_ov_k1 = 2.58e-5;
L.pbranch_ov_k2 = 2.219;
L.pbranch_nov_k1 = 2.34e-3;
L.pbranch_nov_k2 = 0.194;
// Absolute maximum dendritic extent [µm].
L.max_extent = 2000;
return L;
}
lsys_param alpha_motoneuron_lsys = make_alpha_motoneuron_lsys();
// Some parameters missing from the literature are taken from
// the `purkburk.prm` L-Neuron parameter file:
// http://krasnow1.gmu.edu/cn3/L-Neuron/database/purk/bu/purkburk.prm
// Note 'Pov' and 'Pterm' numbers are swapped in Ascoli fig 3B.
lsys_param make_purkinje_lsys() {
lsys_param L;
// Soma diameter [µm].
// (Gaussian fit to 3 samples rom Rapp 1994); truncate to 2σ.
L.diam_soma = { 22.12, 27.68, 24.9, 1.39 };
// Number of dendritic trees (rounded to nearest int).
L.n_tree = { 1 };
// Initial dendrite diameter [µm]. (Dstem)
L.diam_initial = { 4.8, 7.6, 6.167, 1.069 };
// Dendrite step length [µm]. (ΔL)
L.length_step = { 2.3 }; // from `purkburk.prm`
// Tortuousness: roll within section over ΔL [degrees].
L.roll_section = { 0 };
// Tortuousness: pitch within section over ΔL [degrees].
L.pitch_section = { -5, 5 }; // from `purkburk.prm`
// Taper rate: diameter decrease per unit length.
L.taper = { -inf, inf, -0.010, 0.022 };
// Branching torsion: roll at branch point [degrees].
L.roll_at_branch = { -inf, inf, 0.0, 6.5 }; // from `purkburj.prm`
// Branch angle between siblings [degrees].
L.branch_angle = { 1.0, 179.0, 71.0, 40.0 };
// Correlated child branch radius distribution parameters.
L.diam_child_a = 5.47078; // from `purkburk.prm`
L.diam_child_r = { 0, inf, 0.112, 0.035 };
// Termination probability parameters [1/µm].
L.pterm_k1 = 0.1973; // from `purkburj.prm`; 0.4539 in Ascoli.
L.pterm_k2 = -1.1533;
// Bifurcation probability parameters [1/μm].
L.pbranch_ov_k1 = 0.1021; // from `purkburk.prm`; 0.0355 in Ascoli.
L.pbranch_ov_k2 = 0.7237;
L.pbranch_nov_k1 = 0.1021; // from `purkburk.prm`; 0.2349 in Ascoli.
L.pbranch_nov_k2 = -0.0116;
// Absolute maximum dendritic extent [µm].
L.max_extent = 20000;
return L;
}
lsys_param purkinje_lsys = make_purkinje_lsys();
lsys_param make_apical_lsys() {
lsys_param L;
// Soma diameter [µm].
L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
// Number of dendritic trees (rounded to nearest int).
L.n_tree = { 1 }; // numtree
// Initial dendrite diameter [µm]. (Dstem)
L.diam_initial = { 10, 20 }; // initdiam
// Dendrite step length [µm]. (ΔL)
L.length_step = { 25 };
// Initial roll (intrinsic rotation about x-axis) [degrees].
L.roll_initial = { 0 };
// Initial pitch (intrinsic rotation about y-axis) [degrees].
L.pitch_initial = { 0 }; // biforient
// Tortuousness: roll within section over ΔL [degrees].
L.roll_section = { -180.0, 180.0 };
// Tortuousness: pitch within section over ΔL [degrees].
L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
// Taper rate: diameter decrease per unit length.
L.taper = { -1.25e-3 };
// Branching torsion: roll at branch point [degrees].
L.roll_at_branch = { 0.0, 180.0 };
// Branch angle between siblings [degrees].
L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
// Correlated child branch radius distribution parameters.
L.diam_child_a = -0.2087;
L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
// Termination probability parameters [1/µm].
L.pterm_k1 = 2.62e-2;
L.pterm_k2 = -2.955;
// Bifurcation probability parameters [1/μm].
L.pbranch_ov_k1 = 0.58e-5;
L.pbranch_ov_k2 = 2.219;
L.pbranch_nov_k1 = 2.34e-3;
L.pbranch_nov_k2 = 0.194;
// Absolute maximum dendritic extent [µm].
L.max_extent = 2000;
L.tag = 4;
return L;
}
lsys_param apical_lsys = make_apical_lsys();
lsys_param make_basal_lsys() {
lsys_param L;
// Soma diameter [µm].
L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
// Number of dendritic trees (rounded to nearest int).
L.n_tree = { 1, 20, 6, 2}; // numtree
// Initial dendrite diameter [µm]. (Dstem)
L.diam_initial = { 3.0, 12.0 }; // initdiam
// Dendrite step length [µm]. (ΔL)
L.length_step = { 25 };
// Initial roll (intrinsic rotation about x-axis) [degrees].
L.roll_initial = { -180, 180 };
// Initial pitch (intrinsic rotation about y-axis) [degrees].
L.pitch_initial = { 140, 220 }; // biforient
// Tortuousness: roll within section over ΔL [degrees].
L.roll_section = { -180.0, 180.0 };
// Tortuousness: pitch within section over ΔL [degrees].
L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
// Taper rate: diameter decrease per unit length.
L.taper = { -1.25e-3 };
// Branching torsion: roll at branch point [degrees].
L.roll_at_branch = { 0.0, 180.0 };
// Branch angle between siblings [degrees].
L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
// Correlated child branch radius distribution parameters.
L.diam_child_a = -0.2087;
L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
// Termination probability parameters [1/µm].
L.pterm_k1 = 2.62e-2;
L.pterm_k2 = -2.955;
// Bifurcation probability parameters [1/μm].
L.pbranch_ov_k1 = 1.58e-5;
L.pbranch_ov_k2 = 2.219;
L.pbranch_nov_k1 = 2.34e-3;
L.pbranch_nov_k2 = 0.194;
// Absolute maximum dendritic extent [µm].
L.max_extent = 1000;
L.tag = 3;
return L;
}
lsys_param basal_lsys = make_basal_lsys();
#pragma once
#include "lsystem.hpp"
// Predefined parameters for two classes of neurons.
// See lsys_models.cc for details.
extern lsys_param alpha_motoneuron_lsys;
extern lsys_param apical_lsys;
extern lsys_param basal_lsys;
extern lsys_param purkinje_lsys;
#include <cmath>
#include <iostream>
#include <random>
#include <stack>
#include <vector>
#include <map>
#include <arbor/math.hpp>
#include "lsystem.hpp"
using namespace arb::math;
// L-system implementation.
// Random distribution used in the specification of L-system parameters.
// It can be a constant, uniform over an interval, or truncated normal.
// Note that mixture distributions of the above should be, but isn't yet,
// implemented.
struct lsys_distribution {
using result_type = double;
using param_type = lsys_distribution_param;
lsys_distribution(const param_type& p): p_(p) { init(); }
lsys_distribution(): p_() { init(); }
lsys_distribution(result_type x): p_(x) { init(); }
lsys_distribution(result_type lb, result_type ub): p_(lb, ub) { init(); }
lsys_distribution(result_type lb, result_type ub, result_type mu, result_type sigma):
p_(lb, ub, mu, sigma)
{
init();
}
param_type param() const { return p_; }
void param(const param_type& pbis) { p_=pbis; }
result_type min() const { return p_.lb; }
result_type max() const {
return p_.kind==param_type::constant? min(): p_.ub;
}
bool operator==(const lsys_distribution& y) const {
return p_==y.p_;
}
bool operator!=(const lsys_distribution& y) const {
return p_!=y.p_;
}
// omitting out ostream, istream methods
// (...)
template <typename Gen>
result_type operator()(Gen& g, const param_type &p) const {
lsys_distribution dist(p);
return dist(g);
}
template <typename Gen>
result_type operator()(Gen& g) const {
switch (p_.kind) {
case param_type::constant:
return p_.lb;
case param_type::uniform:
return uniform_dist(g);
case param_type::normal:
// TODO: replace with a robust truncated normal generator!
for (;;) {
result_type r = normal_dist(g);
if (r>=p_.lb && r<=p_.ub) return r;
}
}
return 0;
}
private:
param_type p_;
mutable std::uniform_real_distribution<double> uniform_dist;
mutable std::normal_distribution<double> normal_dist;
void init() {
switch (p_.kind) {
case param_type::uniform:
uniform_dist = std::uniform_real_distribution<double>(p_.lb, p_.ub);
break;
case param_type::normal:
normal_dist = std::normal_distribution<result_type>(p_.mu, p_.sigma);
break;
default: ;
}
}
};
class burke_correlated_r {
double a_;
lsys_distribution r_;
public:
struct result_type { double rho1, rho2; };
burke_correlated_r(double a, lsys_distribution_param r): a_(a), r_(r) {}
template <typename Gen>
result_type operator()(Gen& g) const {
double r1 = r_(g);
double r2 = r_(g);
return { r1+a_*r2, r2+a_*r1 };
}
};
class burke_exp_test {
double k1_;
double k2_;
public:
using result_type = int;
burke_exp_test(double k1, double k2): k1_(k1), k2_(k2) {}
template <typename Gen>
result_type operator()(double delta_l, double radius, Gen& g) const {
std::bernoulli_distribution B(delta_l*k1_*std::exp(k2_*radius*2.0));
return B(g);
}
};
class burke_exp2_test {
double ak1_;
double ak2_;
double bk1_;
double bk2_;
public:
using result_type = int;
burke_exp2_test(double ak1, double ak2, double bk1, double bk2):
ak1_(ak1), ak2_(ak2), bk1_(bk1), bk2_(bk2) {}
template <typename Gen>
result_type operator()(double delta_l, double radius, Gen& g) const {
double diam = 2.*radius;
std::bernoulli_distribution B(
delta_l*std::min(ak1_*std::exp(ak2_*diam), bk1_*std::exp(bk2_*diam)));
return B(g);
}
};
// L-system parameter set with instantiated distributions.
struct lsys_sampler {
// Create distributions once from supplied parameter set.
explicit lsys_sampler(const lsys_param& P):
diam_soma(P.diam_soma),
n_tree(P.n_tree),
diam_initial(P.diam_initial),
length_step(P.length_step),
roll_initial(P.roll_initial),
pitch_initial(P.pitch_initial),
roll_section(P.roll_section),
pitch_section(P.pitch_section),
taper(P.taper),
roll_at_branch(P.roll_at_branch),
branch_angle(P.branch_angle),
diam_child(P.diam_child_a, P.diam_child_r),
term_test(P.pterm_k1, P.pterm_k2),
branch_test(P.pbranch_ov_k1, P.pbranch_ov_k2,
P.pbranch_nov_k1, P.pbranch_nov_k2),
max_extent(P.max_extent)
{}
lsys_distribution diam_soma;
lsys_distribution n_tree;
lsys_distribution diam_initial;
lsys_distribution length_step;
lsys_distribution roll_initial;
lsys_distribution pitch_initial;
lsys_distribution roll_section;
lsys_distribution pitch_section;
lsys_distribution taper;
lsys_distribution roll_at_branch;
lsys_distribution branch_angle;
burke_correlated_r diam_child;
burke_exp_test term_test;
burke_exp2_test branch_test;
double max_extent;
};
struct section_point {
double x;
double y;
double z;
double r;
friend std::ostream& operator<<(std::ostream& os, const section_point& obj);
};
std::ostream& operator<<(std::ostream& os, const section_point& obj) {
os << obj.x << " " << obj.y << " " << obj.z << " " << obj.r << '\n';
return os;
}
struct section_tip {
section_point point = {0., 0., 0., 0.};
quaternion rotation = {1., 0., 0., 0.};
double somatic_distance = 0.;
};
struct section_geometry {
unsigned id = 0; // ids should be contigously numbered from 1 in the morphology.
unsigned parent_id = 0;
bool terminal = false;
std::vector<section_point> points;
double length = 0; // µm
section_geometry() = default;
section_geometry(unsigned id, unsigned parent_id, bool terminal, std::vector<section_point> points, double length) :
id(id), parent_id(parent_id), terminal(terminal), points(std::move(points)), length(length)
{}
friend std::ostream& operator<<(std::ostream& os, const section_geometry& obj);
};
std::ostream& operator<<(std::ostream& os, const section_geometry& obj) {
os << "id: " << obj.id << '\n';
os << "parent_id: " << obj.parent_id << '\n';
os << "terminal: " << obj.terminal << '\n';
os << "length: " << obj.length << '\n';
os << "points: " << '\n';
for (auto& p : obj.points) {
os << p << '\n';
}
return os;
}
struct grow_result {
std::vector<section_point> points;
std::vector<section_tip> children;
double length = 0.;
};
constexpr double deg_to_rad = 3.1415926535897932384626433832795/180.0;
template <typename Gen>
grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
constexpr quaternion xaxis = {0, 1, 0, 0};
static std::uniform_real_distribution<double> U;
grow_result result;
std::vector<section_point>& points = result.points;
points.push_back(tip.point);
for (;;) {
quaternion step = xaxis^tip.rotation;
double dl = S.length_step(g);
tip.point.x += dl*step.x;
tip.point.y += dl*step.y;
tip.point.z += dl*step.z;
tip.point.r += dl*0.5*S.taper(g);
tip.somatic_distance += dl;
result.length += dl;
if (tip.point.r<0) tip.point.r = 0;
double phi = S.roll_section(g);
double theta = S.pitch_section(g);
tip.rotation *= rotation_x(deg_to_rad*phi);
tip.rotation *= rotation_y(deg_to_rad*theta);
points.push_back(tip.point);
if (tip.point.r==0 || tip.somatic_distance>=S.max_extent) {
return result;
}
if (S.branch_test(dl, tip.point.r, g)) {
auto r = S.diam_child(g);
// ignore branch if we get non-positive radii
if (r.rho1>0 && r.rho2>0) {
double branch_phi = S.roll_at_branch(g);
double branch_angle = S.branch_angle(g);
double branch_theta1 = U(g)*branch_angle;
double branch_theta2 = branch_theta1-branch_angle;
tip.rotation *= rotation_x(deg_to_rad*branch_phi);
section_tip t1 = tip;
t1.point.r = t1.point.r * r.rho1;
t1.rotation *= rotation_y(deg_to_rad*branch_theta1);
section_tip t2 = tip;
t2.point.r = t2.point.r * r.rho2;
t2.rotation *= rotation_y(deg_to_rad*branch_theta2);
result.children = {t1, t2};
return result;
}
}
if (S.term_test(dl, tip.point.r, g)) {
return result;
}
}
}
std::vector<section_geometry> generate_sections(double soma_radius, lsys_param P, lsys_generator &g) {
constexpr quaternion xaxis = {0, 1, 0, 0};
lsys_sampler S(P);
struct section_start {
section_tip tip;
unsigned parent_id;
};
std::stack<section_start> starts;
unsigned next_id = 1u;
int n = (int)std::round(S.n_tree(g));
for (int i=0; i<n; ++i) {
double phi = S.roll_initial(g);
double theta = S.pitch_initial(g);
double radius = 0.5*S.diam_initial(g);
section_tip tip;
tip.rotation = rotation_x(deg_to_rad*phi)*rotation_y(deg_to_rad*theta);
tip.somatic_distance = 0.0;
auto p = (soma_radius*xaxis)^tip.rotation;
tip.point = {p.x, p.y, p.z, radius};
starts.push({tip, 0u});
}
std::vector<section_geometry> sections;
while (!starts.empty() && sections.size() < P.max_sections) {
auto start = starts.top();
starts.pop();
auto branch = grow(start.tip, S, g);
section_geometry section{next_id++, start.parent_id, branch.children.empty(), std::move(branch.points), branch.length};
for (auto child: branch.children) {
starts.push({child, section.id});
}
sections.push_back(std::move(section));
}
return sections;
}
arb::segment_tree generate_morphology(const lsys_distribution_param& soma, std::vector<lsys_param> Ps, lsys_generator &g) {
double const soma_diameter = lsys_distribution(soma)(g);
double const soma_radius = 0.5*soma_diameter;
arb::segment_tree morph;
morph.append(
arb::mnpos, arb::mpoint{-soma_diameter, 0, 0, soma_diameter},
arb::mpoint{soma_diameter, 0, 0, soma_diameter}, 1);
for(auto P : Ps) {
auto sections = generate_sections(soma_radius, P, g);
// map from internal representation to Arbor ids
std::map<size_t, arb::msize_t> parent_end_id;
// soma
parent_end_id[0] = 0;
// dendrites:
for (auto& sec: sections) {
// the first section must have the soma as parent
size_t parent = parent_end_id.at(sec.parent_id);
const auto& points = sec.points;
if (points.size() < 2) {
throw std::runtime_error("segment has only 1 point");
}
// Include first point only for dendrites segments attached to soma.
if (sec.parent_id==0) {
parent = morph.append(parent,
arb::mpoint{points[0].x, points[0].y, points[0].z, points[0].r},
arb::mpoint{points[1].x, points[1].y, points[1].z, points[1].r}, P.tag);
}
// Remaining points.
for (unsigned i = 1; i < points.size(); ++i) {
const auto& p = points[i];
parent = morph.append(parent, arb::mpoint{p.x, p.y, p.z, p.r}, P.tag);
}
parent_end_id[sec.id] = parent;
}
}
return morph;
}
#pragma once
#include <random>
#include <arbor/morph/segment_tree.hpp>
struct lsys_param;
using lsys_generator = std::minstd_rand;
class lsys_distribution_param;
arb::segment_tree generate_morphology(const lsys_distribution_param& soma, std::vector<lsys_param> Ps, lsys_generator& g);
// The distribution parameters used in the specification of the L-system parameters.
// The distribution can be a constant, uniform over an interval, or truncated normal.
// Note that mixture distributions of the above should be, but isn't yet,
// implemented.
struct lsys_distribution_param {
using result_type = double;
enum lsys_kind {
constant, uniform, normal
};
lsys_kind kind;
result_type lb; // lower bound
result_type ub; // upper bound (used for only non-constant)
result_type mu; // mean (used only for truncated Gaussian)
result_type sigma; // s.d. (used only for truncated Gaussian)
bool operator==(const lsys_distribution_param &p) const {
if (kind!=p.kind) return false;
switch (kind) {
case constant:
return lb == p.lb;
case uniform:
return lb == p.lb && ub == p.ub;
case normal:
return lb == p.lb && ub == p.ub && mu == p.mu && sigma == p.sigma;
}
return false;
}
bool operator!=(const lsys_distribution_param &p) const {
return !(*this==p);
}
// zero or one parameter => constant
lsys_distribution_param():
kind(constant), lb(0.0) {}
lsys_distribution_param(result_type x):
kind(constant), lb(x) {}
// two paramter => uniform
lsys_distribution_param(result_type lb, result_type ub):
kind(uniform), lb(lb), ub(ub) {}
// four paramter => truncated normal
lsys_distribution_param(result_type lb, result_type ub, result_type mu, result_type sigma):
kind(normal), lb(lb), ub(ub), mu(mu), sigma(sigma) {}
};
// Parameters as referenced in Ascoli 2001 are given in parentheses below.
// Defaults below are chosen to be safe, boring and not biological.
struct lsys_param {
// Soma diameter [µm].
lsys_distribution_param diam_soma = { 20 };
// Number of dendritic trees (rounded to nearest int). (Ntree)
lsys_distribution_param n_tree = { 1 };
// Initial dendrite diameter [µm]. (Dstem)
lsys_distribution_param diam_initial = { 2 };
// Dendrite step length [µm]. (ΔL)
lsys_distribution_param length_step = { 25 };
// Dendrites grow along x-axis after application of (instrinsic) rotations;
// roll (about x-axis) is applied first, followed by pitch (about y-axis).
// Initial roll (intrinsic rotation about x-axis) [degrees]. (Taz)
lsys_distribution_param roll_initial = { -45, 45 };
// Initial pitch (intrinsic rotation about y-axis) [degrees]. (Tel)
lsys_distribution_param pitch_initial = { -45, 45 };
// Tortuousness: roll within section over ΔL [degrees]. (Eaz)
lsys_distribution_param roll_section = { -45, 45 };
// Tortuousness: pitch within section over ΔL [degrees]. (Eel)
lsys_distribution_param pitch_section = { -45, 45 };
// Taper rate: diameter decrease per unit length. (TPRB)
lsys_distribution_param taper = { 0 };
// Branching torsion: roll at branch point [degrees]. (Btor)
lsys_distribution_param roll_at_branch = { 0 };
// Branch angle between siblings [degrees]. (Bamp)
lsys_distribution_param branch_angle = { 45 };
// Child branch diameter ratios are typically anticorrelated.
// The Burke 1992 parameterization describes the two child ratios
// as d1 = r1 + a*r2 and d2 = r2 + a*r1 where a is a constant
// determined by the correlation, and r1 and r2 are drawn from
// a common distribution (`d_child_r` below).
double diam_child_a = 0;
lsys_distribution_param diam_child_r = { 0.5 };
// Bifurcation and termination probabilities per unit-length below are given by
// P = k1 * exp(k2*diameter), described by the paremeters k1 and k2 [1/µm].
// Termination probability parameters [1/µm]. (k1trm, k2trm)
double pterm_k1 = 0.05;
double pterm_k2 = -2;
// Bifurcation probability is given by minimum of two probabilities, `pbranch_ov' and
// `pbranch_nov' below.
double pbranch_ov_k1 = 0.01;
double pbranch_ov_k2 = 0;
double pbranch_nov_k1 = 0.01;
double pbranch_nov_k2 = 0;
// Absolute maximum dendritic extent [µm]. (Forces termination of algorithm)
double max_extent = 2000;
// Absolute maximum number of unbranched sections. (Forces termination of algorithm)
unsigned max_sections = 10000;
size_t tag = 0;
};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment