BioDynaMo  v1.03.58-27764645
neuron_soma.cc
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
16 
17 #include <algorithm>
18 #include <ctime>
19 #include <locale>
20 
21 #include "core/resource_manager.h"
24 #include "neuroscience/param.h"
25 
26 namespace bdm {
27 namespace neuroscience {
28 
29 NeuronSoma::NeuronSoma() = default;
30 
31 NeuronSoma::~NeuronSoma() = default;
32 
33 NeuronSoma::NeuronSoma(const Real3& position) : Base(position) {}
34 
36  Base::Initialize(event);
37 
38  if (event.GetUid() == CellDivisionEvent::kUid) {
39  auto* mother = bdm_static_cast<NeuronSoma*>(event.existing_agent);
40  if (mother->daughters_.size() != 0) {
41  Fatal("NeuronSoma",
42  "Dividing a neuron soma with attached neurites is not supported "
43  "in the default implementation! If you want to change this "
44  "behavior derive from this class and overwrite this method.");
45  }
46  }
47 }
48 
49 void NeuronSoma::Update(const NewAgentEvent& event) {
50  Base::Update(event);
51 
52  if (event.GetUid() == NewNeuriteExtensionEvent::kUid) {
53  const auto& ne_event = static_cast<const NewNeuriteExtensionEvent&>(event);
54  auto* neurite = bdm_static_cast<NeuriteElement*>(event.new_agents[0]);
55 
56  real_t theta = ne_event.theta;
57  real_t phi = ne_event.phi;
58  real_t x_coord = std::sin(theta) * std::cos(phi);
59  real_t y_coord = std::sin(theta) * std::sin(phi);
60  real_t z_coord = std::cos(theta);
61 
62  daughters_.push_back(neurite->GetAgentPtr<NeuriteElement>());
63  daughters_coord_[neurite->GetUid()] = {x_coord, y_coord, z_coord};
64  }
65 
66  // do nothing for CellDivisionEvent or others
67 }
68 
69 void NeuronSoma::CriticalRegion(std::vector<AgentPointer<>>* aptrs) {
70  aptrs->reserve(daughters_.size() + 1);
71  aptrs->push_back(Agent::GetAgentPtr<>());
72  for (auto& daughter : daughters_) {
73  aptrs->push_back(daughter);
74  }
75 }
76 
78  NeuriteElement* prototype) {
79  auto dir = direction + GetPosition();
80  auto angles = Base::TransformCoordinatesGlobalToPolar(dir);
81  auto* param = Simulation::GetActive()->GetParam()->Get<Param>();
82  return ExtendNewNeurite(param->neurite_default_diameter, angles[2], angles[1],
83  prototype);
84 }
85 
87  real_t theta,
88  NeuriteElement* prototype) {
89  if (!prototype) {
90  static NeuriteElement kDefaultNeurite;
91  prototype = &kDefaultNeurite;
92  }
93  NewNeuriteExtensionEvent event(diameter, phi, theta);
94  CreateNewAgents(event, {prototype});
95  return bdm_static_cast<NeuriteElement*>(event.new_agents[0]);
96 }
97 
99  auto it = std::find(std::begin(daughters_), std::end(daughters_), daughter);
100  assert(it != std::end(daughters_) &&
101  "The element you wanted to remove is not part of daughters_");
102  daughters_.erase(it);
103 }
104 
105 Real3 NeuronSoma::OriginOf(const AgentUid& daughter_uid) const {
106  Real3 xyz = daughters_coord_.at(daughter_uid);
107 
108  real_t radius = GetDiameter() * .5;
109  xyz = xyz * radius;
110 
111  Real3 axis_0 = {Base::kXAxis[0], Base::kYAxis[0], Base::kZAxis[0]};
112  Real3 axis_1 = {Base::kXAxis[1], Base::kYAxis[1], Base::kZAxis[1]};
113  Real3 axis_2 = {Base::kXAxis[2], Base::kYAxis[2], Base::kZAxis[2]};
114 
115  Real3 result = {xyz * axis_0, xyz * axis_1, xyz * axis_2};
116  return GetPosition() + result;
117 }
118 
120 
122  const NeuronOrNeurite& new_rel) {
123  auto old_rel_agent_ptr = bdm_static_cast<const NeuriteElement*>(&old_rel)
124  ->GetAgentPtr<NeuriteElement>();
125  auto new_rel_agent_ptr = bdm_static_cast<const NeuriteElement*>(&new_rel)
127  auto coord = daughters_coord_[old_rel_agent_ptr->GetUid()];
128  auto it = std::find(std::begin(daughters_), std::end(daughters_),
129  old_rel_agent_ptr);
130  assert(it != std::end(daughters_) &&
131  "old_element_idx could not be found in daughters_ vector");
132  *it = new_rel_agent_ptr;
133  daughters_coord_[new_rel_agent_ptr->GetUid()] = coord;
134 }
135 
136 const std::vector<AgentPointer<NeuriteElement>>& NeuronSoma::GetDaughters()
137  const {
138  return daughters_;
139 }
140 
141 // Helper function for PrintSWC to print the elements of a vector seperated by
142 // a space.
143 inline void PrintVector(std::ostream& out, const Real3& pos) {
144  out << " " << pos[0] << " " << pos[1] << " " << pos[2] << " ";
145 }
146 
147 void NeuronSoma::PrintSWC(std::ostream& out) const {
148  // 1. Write a header to the file such that users do not confuse it with real
149  // data.
150  out << "# This SWC file was created with BioDynaMo and is (very) likely "
151  "to\n# be the result of a simulation. Be careful not to confuse it\n"
152  "# with real data. \n\n";
153  // Example for localtime_r: https://tinyurl.com/localtimer
154  struct tm newtime;
155  time_t ltime;
156  ltime = time(&ltime);
157  localtime_r(&ltime, &newtime);
158  char creation_time[100];
159  std::strftime(creation_time, sizeof(creation_time), "%A, %F, %T", &newtime);
160  out << "# Created at: " << creation_time << "\n\n";
161 
162  // 2. Write the line for the soma at the beginning of the file
163  int element_id{1};
164  int previous_element_id{-1};
165  out << element_id++ << " " << GetIdentifierSWC();
166  PrintVector(out, GetPosition());
167  out << GetDiameter() / 2 << " " << previous_element_id << "\n";
168 
169  // 3. Iterate over all neuron branches that are attached to the soma
170  for (auto daughter : daughters_) {
171  // New branches need to be connected to soma
172  bool is_new_brach{true};
173  // Pointers to currently considered and next element
174  AgentPointer<NeuriteElement> current_element, next_element;
175  // Vector to track the bifurcation elements (right daughters that are
176  // neglected when going down the left path, see below)
177  std::vector<AgentPointer<NeuriteElement>> bifurcation_points;
178  // Track the labels that we need to connect the bifurcation points to the
179  // appropriate elements
180  std::vector<int> bifurcation_points_lables;
181  // Used to reconnect to previous bifurcation
182  bool connect_to_bifurcation{false};
183 
184  // The logic of the while loop is the following. The tree of dendrites /
185  // axons is made up of subsequent elements. Each element has a left daughter
186  // but only bifurcations points also have a right daughter. If a branch
187  // ends, the element has neither a left nor a right daughter. Thus, for each
188  // branch attached to the soma, we walk down the tree all the way down
189  // following the left path. On the way, we remember the bifurcation points.
190  // Once we've reached the end of a branch, we jump back to the last
191  // bifurcation point that we have visited and start to follow the left path
192  // again. With this method we can resolve the entire tree structure under
193  // the assumption, that each element has 0, 1, or 2 follow up elements and
194  // that different branches do not merge again.
195  current_element = daughter;
196  while (true) {
197  // Loop 1. Write down the SWC line corresponding to the current element.
198  // The file format per line is as follows:
199  // <element id> <type id> <x> <y> <z> <radius> <prev element id>
200  out << element_id << " " << current_element->GetIdentifierSWC();
201  PrintVector(out, current_element->GetPosition());
202  out << current_element->GetDiameter() / 2 << " ";
203  // This block establishes the connection via <prev element id>
204  if (is_new_brach) {
205  // New branches need to be attached to the soma
206  previous_element_id = 1;
207  is_new_brach = false;
208  } else {
209  if (!connect_to_bifurcation) {
210  // By default we connect to the previously considered element
211  previous_element_id = element_id - 1;
212  } else {
213  // At some points we need to reconnect the branches to the previous
214  // bifurcation point.
215  previous_element_id = bifurcation_points_lables.back();
216  bifurcation_points_lables.pop_back();
217  connect_to_bifurcation = false;
218  }
219  }
220  out << previous_element_id << "\n";
221 
222  // Loop 2. Walk down one side of the tree, take left option first, and
223  // remember bifurcation points.
224  next_element = current_element->GetDaughterLeft();
225  // Track the bifurcations
226  if (current_element->GetDaughterRight() != nullptr) {
227  bifurcation_points.push_back(current_element->GetDaughterRight());
228  bifurcation_points_lables.push_back(element_id);
229  }
230  element_id++; // Increase ID by 1
231  if (next_element == nullptr && bifurcation_points.empty()) {
232  // If we have resolved all bifurcations and don't have any next element,
233  // we exit the loop
234  break;
235  } else if (next_element == nullptr) {
236  // Jump to previous bifurcation
237  next_element = bifurcation_points.back();
238  bifurcation_points.pop_back();
239  connect_to_bifurcation = true;
240  }
241  std::swap(current_element, next_element); // Swap Pointers
242  }
243  }
244 }
245 
246 } // namespace neuroscience
247 } // namespace bdm
bdm::NewAgentEvent
Definition: new_agent_event.h:61
bdm::neuroscience::NeuronSoma::daughters_coord_
std::unordered_map< AgentUid, Real3 > daughters_coord_
Definition: neuron_soma.h:128
bdm::neuroscience::NeuronSoma::~NeuronSoma
virtual ~NeuronSoma()
bdm::neuroscience::NewNeuriteExtensionEvent
Contains the parameters to extend a new neurite from a neuron soma.
Definition: new_neurite_extension_event.h:27
bdm::kYAxis
@ kYAxis
Definition: substance_initializers.h:33
bdm::neuroscience::PrintVector
void PrintVector(std::ostream &out, const Real3 &pos)
Definition: neuron_soma.cc:143
bdm
Definition: agent.cc:39
bdm::neuroscience::NeuronSoma::Update
void Update(const NewAgentEvent &event) override
This method is used to update attributes after a cell division. or new neurite branching event.
Definition: neuron_soma.cc:49
bdm::Agent::CreateNewAgents
void CreateNewAgents(const NewAgentEvent &event, const std::initializer_list< Agent * > &prototypes)
Definition: agent.h:121
bdm::AgentPointer
Definition: agent_pointer.h:58
bdm::neuroscience::NeuronSoma::Initialize
void Initialize(const NewAgentEvent &event) override
This method is used to initialise the values of daughter 2 for a cell division event.
Definition: neuron_soma.cc:35
bdm::real_t
double real_t
Definition: real_t.h:21
neuron_soma.h
bdm::neuroscience::NeuronSoma::NeuronSoma
NeuronSoma()
new_neurite_extension_event.h
bdm::kZAxis
@ kZAxis
Definition: substance_initializers.h:33
neurite_element.h
bdm::Cell::GetPosition
const Real3 & GetPosition() const override
Definition: cell.h:188
bdm::neuroscience::NeuronSoma::PrintSWC
void PrintSWC(std::ostream &out) const
Exports the soma and the attached neurite elements to the SWC file format.
Definition: neuron_soma.cc:147
bdm::NewAgentEvent::existing_agent
Agent * existing_agent
Definition: new_agent_event.h:69
bdm::Param::Get
const TParamGroup * Get() const
Definition: param.h:59
bdm::neuroscience::Param
Definition: param.h:32
bdm::neuroscience::NeuronSoma::UpdateDependentPhysicalVariables
void UpdateDependentPhysicalVariables() override
Definition: neuron_soma.cc:119
bdm::neuroscience::NeuriteElement::GetUid
const AgentUid & GetUid() const override
Definition: neurite_element.h:66
bdm::neuroscience::NeuronOrNeurite
Definition: neuron_or_neurite.h:42
bdm::neuroscience::NeuronOrNeurite::GetIdentifierSWC
virtual StructureIdentifierSWC GetIdentifierSWC() const
Returns the SWC classification of the object.
Definition: neuron_or_neurite.cc:35
bdm::NewAgentEvent::new_agents
InlineVector< Agent *, 3 > new_agents
Definition: new_agent_event.h:73
param.h
bdm::neuroscience::NeuronSoma::daughters_
std::vector< AgentPointer< NeuriteElement > > daughters_
Definition: neuron_soma.h:123
bdm::neuroscience::NeuronSoma::CriticalRegion
void CriticalRegion(std::vector< AgentPointer<>> *aptrs) override
Definition: neuron_soma.cc:69
bdm::neuroscience::NeuronSoma::UpdateRelative
void UpdateRelative(const NeuronOrNeurite &old_rel, const NeuronOrNeurite &new_rel) override
Definition: neuron_soma.cc:121
bdm::Cell::GetDiameter
real_t GetDiameter() const override
Definition: cell.h:182
bdm::Agent::GetAgentPtr
AgentPointer< TAgent > GetAgentPtr() const
Return agent pointer.
Definition: agent.h:208
bdm::AgentUid
Definition: agent_uid.h:25
bdm::neuroscience::NeuronSoma::GetDaughters
const std::vector< AgentPointer< NeuriteElement > > & GetDaughters() const
Definition: neuron_soma.cc:136
bdm::NewAgentEvent::GetUid
virtual NewAgentEventUid GetUid() const =0
bdm::neuroscience::NeuronSoma::OriginOf
Real3 OriginOf(const AgentUid &daughter_uid) const override
Definition: neuron_soma.cc:105
bdm::Simulation::GetParam
const Param * GetParam() const
Returns the simulation parameters.
Definition: simulation.cc:271
bdm::MathArray
Definition: math_array.h:36
resource_manager.h
bdm::Simulation::GetActive
static Simulation * GetActive()
This function returns the currently active Simulation simulation.
Definition: simulation.cc:67
bdm::CellDivisionEvent::kUid
static const NewAgentEventUid kUid
Definition: cell_division_event.h:32
bdm::neuroscience::NewNeuriteExtensionEvent::kUid
static const NewAgentEventUid kUid
Definition: new_neurite_extension_event.h:28
bdm::neuroscience::NeuronSoma::RemoveDaughter
void RemoveDaughter(const AgentPointer< NeuriteElement > &daughter) override
Definition: neuron_soma.cc:98
bdm::kXAxis
@ kXAxis
Definition: substance_initializers.h:33
bdm::neuroscience::NeuronSoma::ExtendNewNeurite
NeuriteElement * ExtendNewNeurite(const Real3 &direction, NeuriteElement *prototype=nullptr)
Extend a new neurite from this soma.
Definition: neuron_soma.cc:77
bdm::neuroscience::NeuriteElement
Definition: neurite_element.h:56