Моделирование электричесой активности нейронов C# реализация

Данный класс представляет из себя реализацию модели описанной в статье автора esir_pavel

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
 
namespace SpikeNN
{
    class SpikeNN
    {
        const float h     = .5f;    // временной шаг интегрирования в мс
        const int   Tsim  = 500; // время симуляции в дискретных отсчетах
        const int   Nexc  = 100;    // Количество возбуждающих (excitatory) нейронов
        const int   Ninh  = 25;     // Количество тормозных (inhibitory) нейронов
        const int   Nneur = Nexc + Ninh;
        const int   Ncon  = (int)(Nneur*Nneur*0.1f); // Количество сязей, 0.1 это вероятность связи между 2-мя случайными нейронами
 
        float[,] Vms = new float[Nneur,Tsim]; // мембранные потенциалы
        float[,] Ums = new float[Nneur,Tsim]; // вспомогательные переменные модели Ижикевича
        float[] Iex = new float [Nneur];   // внешний постоянный ток приложенный к нейрону
        float[] Isyn= new float [Nneur];  // синаптический ток на каждый нейрон
 
        int[] pre_conns= new int [Ncon];     // индексы пресинаптических нейронов
        int[] post_conns= new int [Ncon];   // индексы постсинаптических нейронов
        float[] weights = new float[Ncon];  // веса связей
        float[,] y = new float[Ncon,Tsim];    // переменная модулирующая синаптический ток в зависимости от спайков на пресинапсе
 
        float psc_excxpire_time = 4.0f; // характерное вермя спадания постсинаптического тока, мс
        float minWeight = 50.0f; // веса, размерность пкА
        float maxWeight = 100.0f;
 
        // Параметры нейрона
        float Iex_max = 40.0f; // максимальный приложенный к нейрону ток 50 пкА
        float a    = 0.02f;
        float b    = 0.5f;
        float c    = -40.0f; // значение мембранного потенциала до которого он сбрасываеться после спайка
        float d    = 100.0f;
        float k    = 0.5f;
        float Vr   = -60.0f;
        float Vt   = -45.0f;
        float Vpeak   = 35.0f;  // максимальное значение мембранного потенциала, при котором происходит сброс до значения с
        float V0   = -60.0f; // начальное значение для мембранного потенциала
        float U0   = 0.0f;   // начальное значение для вспомогательной переменной
        float Cm      = 50.0f;  // электрическая ёмкость нейрона, размерность пкФ
 
        float[] spike_times=new float[Nneur*Tsim];    // времена возникновения спайков
        int[] spike_neurons = new int[Nneur * Tsim];  // индексы нейронов на которых происходят спайки
        int spike_num = 0;
 
        Random rnd = new Random();
        void init_connections()
        {
            for (int con_idx = 0; con_idx < Ncon; )
            {
                // случайно выбираем постсипантические и пресинаптические нейроны
                int pre = rnd.Next(Nneur);
                int post = rnd.Next(Nneur);
                pre_conns[con_idx] = pre;
                post_conns[con_idx] = post;
                weights[con_idx] = rnd.Next((int)(((int)(maxWeight - minWeight) * 10) / 10.0f)) + minWeight;
                if (pre >= Nexc)
                {
                    // если пресинаптический нейрон тормозный то вес связи идет со знаком минус
                    weights[con_idx] = -weights[con_idx];
                }
                y[con_idx,0] = 0.0f;
                con_idx++;
            }
        }
 
        void init_neurons()
        {
            for (int neur_idx = 0; neur_idx < Nneur; neur_idx++)
            {
                // случайно разбрасываем приложенные токи
                Iex[neur_idx] = (int)(rnd.Next( (int)(Iex_max * 10)) / 10.0f);
                Isyn[neur_idx] = 0.0f;
                Vms[neur_idx,0] = V0;
                Ums[neur_idx,0] = U0;
            }
        }
 
 
        float izhik_Vm(int neuron, int time)
        {
            return (k * (Vms[neuron,time] - Vr) * (Vms[neuron,time] - Vt) - Ums[neuron,time] + Iex[neuron] + Isyn[neuron]) / Cm;
        }
 
        float izhik_Um(int neuron, int time)
        {
            return a * (b * (Vms[neuron,time] - Vr) - Ums[neuron,time]);
        }
 
        void save2file()
        {
            List<string> rastr = new List<string>();
            string s="";
            for (int k = 0; k < spike_num; k++)
            {
                s=spike_times[k]+"; " + (spike_neurons[k] + 1).ToString() + "; \n" ;
                rastr.Add(s);
            }
            File.WriteAllLines("rastr.csv", rastr);
 
            // Вычисление среднего по всей сети мембранного потенциала в каждый момент времени
            // нечто наподобие электроэнцефалографии
            List<string> oscill = new List<string>();
            for (int t = 0; t < Tsim; t++)
            {
                float Vm_mean = 0.0f;
                for (int m = 0; m < Nneur; m++)
                {
                    Vm_mean += Vms[m,t];
                }
                Vm_mean /= Nneur;
                s = (t * h).ToString() + "; " + Vm_mean.ToString() + "; \n";
                oscill.Add(s);
            }
            File.WriteAllLines("oscill.csv", oscill);
        }
 
       public  void Start()
        {
            init_connections();
            init_neurons();
            float expire_coeff = -h / psc_excxpire_time;
            expire_coeff =(float) Math.Exp(Convert.ToDouble(expire_coeff));
            for (int t = 1; t < Tsim; t++)
            {
                // проходим по всем нейронам
                for (int neur = 0; neur < Nneur; neur++)
                {
                    Vms[neur,t] = Vms[neur,t - 1] + h * izhik_Vm(neur, t - 1);
                    Ums[neur,t] = Ums[neur,t - 1] + h * izhik_Um(neur, t - 1);
                    Isyn[neur] = 0.0f;
 
                    if (Vms[neur,t - 1] > Vpeak)
                    {
                        Vms[neur,t] = c;
                        Ums[neur,t] = Ums[neur,t - 1] + d;
                        spike_times[spike_num] = t * h;
                        spike_neurons[spike_num] = neur;
                        spike_num++;
                    }
                }
 
                // проходим по всем связям
                for (int con = 0; con < Ncon; con++)
                {
                    y[con,t] = y[con,t - 1] * expire_coeff;
 
                    if (Vms[pre_conns[con],t - 1] > Vpeak)
                    {
                        y[con,t] = 1.0f;
                    }
                    Isyn[post_conns[con]] += y[con,t] * weights[con];
                }
            }
            save2file();
        }
    }
}