Выпускная квалификационная работа Процесс испускания вещества в решении задачи диффузии при симуляции растущих нейрональных структур (стр. 4 )

Рисунок 11. Визуализация диффузионного процесса

2.6 Сравнительный анализ

       После всех этапов разработки была произведена симуляция диффузии 100000 узлов с использованием структуры октордерева и массива. В ходе сравнительного анализа была доказана эффективность полученного решения при вычислении больших объемов данных (рисунок 12).

       Из приведенного графика следует, что разработанное программное решение показало лучшие результаты по времени работы.

       

Рисунок 12. Сравнительный анализ производительности алгоритмов. Красная линия - регулярная сетка, синяя линия - структура октодерева

ЗАКЛЮЧЕНИЕ

В результате работы была улучшена существующая структура и разработано программное решение, показывающие лучшие результаты в сравнении с существующими на данный момент аналогами.

При выполнении работы были решены следующие задачи:

Анализ существующих на данный момент решений для симуляции роста нейронных структур; Реализация программного решения для процесса испускания вещества; Сравнение реализованного программного решения с существующим аналогом.

Дальнейшая разработка проекта предполагает разработку алгоритмов октодерева для работы в облаке с использованием параллельных вычислений. Следующим этапом в развитии диффузии веществ будет посвящена реализация существующих алгоритмов с использованием методов параллельных вычислений.

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. The Code for Facial Identity in the Primate Brain. [Текст] – с.1014-1017, 2017 г.

2. Biodynamo: a platform for computer simulations of biological dynamics [Электронный ресурс] - https://biodynamo. web. cern. ch/, режим доступа – свободный.

3. Poplawski NJ, Shirinifard A, Swat M, Glazier JA. Simulation of single-species bacterial-biofilm growth using the Glazier-Graner-Hogeweg model and the CompuCell3D modeling environment. Math Biosci Eng, [Текст]  – с. 13-17, 2008г.

4. F. Zubler and R. Douglas. A framework for modeling the growth and development of neurons and networks. Frontiers in Computational Neuroscience, [Текст]  – с. 10-12, 2009 г.

5. NEURON for empirically-based simulations of neurons and networks of neurons [Электронный ресурс] - http://www. neuron. yale. edu/neuron/, режим доступа – свободный.

6. GENESIS 2 simulator [Электронный ресурс] - http://genesis-sim. org/, режим доступа – свободный.

7. PCSIM: A Parallel neural Circuit SIMulator [Электронный ресурс] - http://www. lsm. tugraz. at/pcsim/ , режим доступа – свободный.

8. Cx3D: Cortex simulation in 3D [Электронный ресурс] - https://www. ini. uzh. ch/~amw/seco/cx3d/, режим доступа – свободный.

9. D. Meagher. Geometric modeling using octree puter Graphics and Image Processing, [Текст]  – c.129–147, 1982 г.

10. Advanced Octrees 2: node representations [Электронный ресурс] - https://geav. /2014/08/18/advanced-octrees-2-node-representations/, режим доступа – свободный.

11. I. Gargantini. An Effective Way to Represent munications of the ACM, [Текст]  – с. 905-910, 1982г.

12. John K. Chilton. Molecular mechanisms of axon guidance. Developmental Biology, [Текст] – c.292 2006 г.

13. leymanov, F. Gafarov, and N. Khusnutdinov. Modeling of interstitial branching of axonal networks. Journal of Integrative Neuroscience, [Текст] – с.- 12, 2013 г.

ПРИЛОЖЕНИЕ

diffusion_op. h

#ifndef BIODYNAMO_DIFFUSION_OP_H

#define BIODYNAMO_DIFFUSION_OP_H

#include <vector>

namespace bdm {

using spatial_organization::SpatialTreeNode;

using spatial_organization::Point;

using spatial_organization::Voxel;

using spatial_organization::AdaptiveMesh;

using std::vector;

using std::array;

// TODO(Sotnem13):

template <typename T>

struct ConcentrationChange {

  Voxel<T>* top;

  Voxel<T>* bottom;

  double delta_value;

};

class DiffusionOp {

public:

  DiffusionOp() {}

  explicit DiffusionOp(int iteration_count = 1, double diffusion_constant = 0.5)

  : iteration_count_(iteration_count), diffusion_constant_(diffusion_constant){}

  ~DiffusionOp() {}

  DiffusionOp(const DiffusionOp&) = delete;

  DiffusionOp& operator=(const DiffusionOp&) = delete;

  template <typename T>

  void Compute(AdaptiveMesh<T>* mesh) const {

  for (int i = 0; i < iteration_count_; i++) {

  mesh->Refine([D = diffusion_constant_]

  (typename AdaptiveMesh<T>::Element &mesh_element) {

  auto neighbors = mesh_element. Neighbors();

  auto delta_value = 0; ;

  for (auto neighbor: neighbors) {

  delta_value += (neighbor. first - mesh_element. Value());

  }

  auto new_value = mesh_element. Value();

  new_value += delta_value*D/6;

  mesh_element. SetValue(new_value);

  });

  }

  }

private:

  int iteration_count_;

  double diffusion_constant_;

};

}  // namespace bdm

#endif //BIODYNAMO_DIFFUSION_OP_H

param. h

#ifndef PARAM_H_

#define PARAM_H_

namespace bdm {

class Param {

public:

  /// Time between two simulation step, in hours.

  static constexpr double kSimulationTimeStep = 0.01;

  /// Maximum jump that a point mass can do in one time step. Useful to

  /// stabilize the simulation

  static constexpr double kSimulationMaximalDisplacement = 3.0;

  static constexpr double kDefaultTolerance = 0.000000001;

  /// Maximum length of a discrete segment before it is cut into two parts.

  static double kNeuriteMaxLength;  // usual value : 20

  /// Minimum length of a discrete segment before. If smaller it will try to

  /// fuse with the proximal one

  static constexpr double kNeuriteMinLength = 2.0;  // usual value : 10

  // Diffusion (saving time by not running diffusion for too small differences)

  /// If concentration of a substance smaller than this, it is not diffused

  static constexpr double kMinimalConcentrationForExtracellularDiffusion = 1e-5;

  /// If absolute concentration difference is smaller than that, there is no

  /// diffusion

  static constexpr double

  kMinimalDifferenceConcentrationForExtracacellularDiffusion = 1e-5;

  /// If ratio (absolute concentration difference)/concentration is smaller than

  /// that, no diffusion.

  static constexpr double kMinimalDCOverCForExtracellularDiffusion = 1e-3;

  /// If concentration of a substance smaller than this, it is not diffused

  static constexpr double kMinimalConcentrationForIntracellularDiffusion =

  1e-10;

  /// If absolute concentration difference is smaller than that, there is no

  /// diffusion

  static constexpr double

  kMinimalDifferenceConcentrationForIntracellularDiffusion = 1e-7;

  /// If ration (absolute concentration difference)/concentration is smaller

  /// than that, no diffusion.

  static constexpr double kMinimalDCOverCForIntracellularDiffusion = 1e-4;

  // Neurites

  /// Initial value of the restingLength before any specification.

  static constexpr double kNeuriteDefaultActualLength = 1.0;

  /// Diameter of an unspecified (= axon/dendrite) neurite when extends from the

  /// somaElement

  static constexpr double kNeuriteDefaultDiameter = 1.0;

  static constexpr double kNeuriteMinimalBifurcationLength = 0;

  /// Spring constant

  static constexpr double kNeuriteDefaultSpringConstant = 10;  // 10;

  /// Threshold the force acting on a neurite has to reach before a move is made

  /// ( = static friction).

  static constexpr double kNeuriteDefaultAdherence = 0.1;

  /// Rest to the movement ( = kinetic friction).

  static constexpr double kNeuriteDefaultMass = 1;

  static constexpr double kNeuriteDefaultTension = 0.0;

  // Somata

  /// CAUTION: not the radius but the diameter

  static constexpr double kSphereDefaultDiameter = 20;

  /// Threshold the force acting on a somaElement has to reach before a move is

  /// made ( = static friction).

  static constexpr double kSphereDefaultAdherence = 0.4;

  /// Restistance to the movement ( = kinetic friction).

  static constexpr double kSphereDefaultMass = 1;

  static constexpr double kSphereDefaultRotationalInertia = 0.5;

  static constexpr double kSphereDefaultInterObjectCoefficient = 0.15;

  /// Helpful constant to compare with 0

  static constexpr double kEpsilon = 1e-20;

  /// Helpful constant to identify 'infinity'

  static constexpr double kInfinity = 1e20;

};

}  // namespace bdm

#endif  // PARAM_H_

dividing_cell_op. h

#ifndef DIVIDING_CELL_OP_H_

#define DIVIDING_CELL_OP_H_

#include "backend. h"

namespace bdm {

class DividingCellOp {

public:

  DividingCellOp() {}

  ~DividingCellOp() {}

  DividingCellOp(const DividingCellOp&) = delete;

  DividingCellOp& operator=(const DividingCellOp&) = delete;

  template <typename TContainer>

  void Compute(TContainer* cells) const {

#pragma omp parallel for

  for (size_t i = 0; i < cells->size(); i++) {

  // if diameter <= 40 then changeVolume(300) else do nothing

  auto&& cell = (*cells)[i];

  if (cell. GetDiameter() <= 40) {

  cell. ChangeVolume(300);

  }

  // todo(lukas) division if diameter > 20;

  }

  }

};

}  // namespace bdm

#endif  // DIVIDING_CELL_OP_H_

cell. h

#ifndef CELL_H_

#define CELL_H_

#include <array>

#include <cmath>

#include <type_traits>

#include "backend. h"

#include "cell. h"

#include "default_force. h"

#include "inline_vector. h"

#include "math_util. h"

#include "param. h"

#include "simulation_object. h"

#include "simulation_object_util. h"

namespace bdm {

using std::array;

template <typename Base = SimulationObject<>>

class CellExt : public Base {

  BDM_CLASS_HEADER(CellExt, position_, mass_location_, tractor_force_,

  diameter_, volume_, adherence_, mass_, neighbors_)

public:

  CellExt() {}

  explicit CellExt(double diameter) : diameter_(diameter) { UpdateVolume(); }

  explicit CellExt(const array<double, 3>& position)

  : position_(position), mass_location_(position) {}

  virtual ~CellExt() {}

  double GetAdherence() const { return adherence_[idx_]; }

  double GetDiameter() const { return diameter_[idx_]; }

  double GetMass() const { return mass_[idx_]; }

  const array<double, 3>& GetMassLocation() const {

  return mass_location_[idx_];

  }

  const array<double, 3>& GetPosition() const { return position_[idx_]; }

  const array<double, 3>& GetTractorForce() const {

  return tractor_force_[idx_];

  }

  double GetVolume() const { return volume_[idx_]; }

  const InlineVector<int, 8>& GetNeighbors() const { return neighbors_[idx_]; }

  void SetAdherence(double adherence) { adherence_[idx_] = adherence; }

  void SetDiameter(double diameter) { diameter_[idx_] = diameter; }

  void SetMass(double mass) { mass_[idx_] = mass; }

  void SetMassLocation(const array<double, 3>& mass_location) {

  mass_location_[idx_] = mass_location;

  }

  void SetPosition(const array<double, 3>& position) {

  position_[idx_] = position;

  }

  void SetTractorForce(const array<double, 3>& tractor_force) {

  tractor_force_[idx_] = tractor_force;

  }

  void SetNeighbors(const InlineVector<int, 8>& neighbors) {

  neighbors_[idx_] = neighbors;

  }

  void ChangeVolume(double speed) {

  // scaling for integration step

  double dV = speed * Param::kSimulationTimeStep;

  volume_[idx_] += dV;

  if (volume_[idx_] < 5.2359877E-7) {

  volume_[idx_] = 5.2359877E-7;

  }

  UpdateDiameter();

  }

  void UpdateDiameter() {

  // V = (4/3)*pi*r^3 = (pi/6)*diameter^3

  diameter_[idx_] = std::cbrt(volume_[idx_] * 6 / Math::kPi);

  }

  void UpdateVolume() {

  volume_[idx_] =

  Math::kPi / 6 * diameter_[idx_] * diameter_[idx_] * diameter_[idx_];

  }

  void UpdateMassLocation(const array<double, 3>& delta) {

  mass_location_[idx_][0] += delta[0];

  mass_location_[idx_][1] += delta[1];

  mass_location_[idx_][2] += delta[2];

  }

  void GetForceOn(const array<double, 3>& ref_mass_location,

  double ref_diameter, array<double, 3>* force) const {

  DefaultForce default_force;

  double iof_coefficient = Param::kSphereDefaultInterObjectCoefficient;

  default_force. ForceBetweenSpheres(ref_mass_location, ref_diameter,

  iof_coefficient, mass_location_[idx_],

  diameter_[idx_], iof_coefficient, force);

  }

private:

  vec<array<double, 3>> position_;

  vec<array<double, 3>> mass_location_;

  vec<array<double, 3>> tractor_force_;

  vec<double> diameter_;

  vec<double> volume_;

  vec<double> adherence_;

  vec<double> mass_;

  // stores a list of neighbor ids for each scalar cell

  vec<InlineVector<int, 8>> neighbors_;

};

template <typename Backend = Scalar>

using Cell = CellExt<SimulationObject<Backend>>;

}  // namespace bdm

#endif  // CELL_H_