Однією з важливих галузей застосування інструментальних систем в медицині залишається медична наука, особливо такий її розділ як математичне моделювання медичних процесів. Більшість з них носить складний нелінійний характер.

Для представлення закономірностей розвитку медичних процесів слід використовувати системи з післядією. Динамічні системи з післядією відіграють важливу роль в економіці, біології, медицині, техніці. Чисельні застосування ередитарних систем показані в роботах [1-3]. У той же час широке практичне використання вказаного класу систем пов’язане з певними труднощами, що виникають при чисельному інтегруванні (розв”язуванні) відповідних диференціальних рівнянь. Однією з причин є, мабуть, нескінченновимірність задач. Нескінченновимірність задачі означає, що не існує скінченного числа деяких базових розв”язків, через які буде виражатися будь-який інший розв”язок системи. У цьому, можливо, і причина недостатнього розвитку аналітичної теорії (теорії явних розв”язків) вказаного класу систем. Тому слід вдатися до чисельно-комп”ютерного розв’язування.

Розглядається клас динамічних систем, що описуються диференціальними рівняннями із запізненням:

  ,     .                      (1)

Тут    – незалежна змінна,  – фазова траєкторія системи,  – змінні запізнення в системі,   – функціонал, який може бути навіть кусково-неперервним.

Припускається, що функції   приймають на додатній півосі невід”ємні значення і обмежені знизу і зверху. Поведінка системи на початковому проміжку часу вважається заданою:

                                                               (2)

Тут  – функція, що допускає розв”язок початкової задачі (1), (2).

Розглянемо метод, придатний для комп’ютерного розв”язування початкових задач загального вигляду (1), (2).

При чисельному інтегруванні звичайних диференціальних рівнянь найбільший розвиток отримав метод Рунге-Кутти і численні його модифікації. Стосовно  задач (1), (2) метод Рунге-Кутти можна застосувати лише у випадку одного сталого запізнення  .

Причому це може бути лише метод з постійною довжиною кроку, що приводить до значних похибок і відзначається поганою стійкістю. При розгляді більш складних систем і використанні більш досконалих методів із змінною довжиною кроку виникають наступні труднощі:

n     при розгляді декількох змінних запізнень не завжди відомо наскільки великий діапазон попередніх станів системи потрібно зберігати для обчислення подальших;

n     не завжди зберігаються саме ті значення, які будуть вимагатися на майбутніх кроках;

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

Всі вищеперелічені труднощі дозволяє уникнути глобальна апроксимація розв”язку. Відповідними методами для щільної видачі розв”язків є методи типу Адамса, неперервні методи Рунге-Кутти. У зв’язку з легкістю програмування і переваги над іншими інтерполяційний методами [3] в даній роботі вибраний метод Дормана і Прінса порядку 5(4).

            Програмна побудова розв”язку задачі (1), (2) будується на використанні ООП в середовищі Java.

            На малюнку 1 наведена схема класів, включених в проект (Workspace) для розв”язування задач (1), (2). Тут DelaySystemSolution – клас, в якому описується початкова поведінка і праві частини системи (1) а також здійснюється чисельне інтегрування методом Дормана-Принса; GraphConstruction – класс-аплет, що здійснює візуальне відображення графіків розв”язків, їх ретельніше відображення і вивчення; FunctionList  – клас – “пов”язаний список”, призначений для розв”язування задачі одночасного виведення кількох компонент розв”язку в одній площині; ColorList – клас, призначений для “різнокольорового” відображення графіків різних компонент розв”язку; LinkedList – клас – “пов”язаний список”, призначений для відображення розв”язку задачі (1), (2) в різних межах зміни  ; BoundsLocation – клас для збереження меж зміни  .

Стрілками позначені напрямки взаємодії класів. Клас GraphConstruction є типовою бібліотекою машинної графіки, що має багаті візуальні можливості щодо вивчення поведінки системи на найдрібніших інтервалах зміни  , одночасного виведення кількох графіків в одній площині.

Клас використовує багатопоточність, а з метою уникнення мерехтіння використовується метод подвійної буферизації. Нижче вказані головні змінні і методи класу GraphConstruction.

public class GraphConstruction extends java.applet.Applet implements Runnable

{

            DelaySystemSolution SampleSystemSolution;             

            Thread m_Graphics = null; // ïîò³ê êëàñó

            private double m_x0;   // ë³âà ìåæà õ

            private double m_x1;   // ïðàâà ìåæà õ

            // другий буфер зображення

            private Image m_image; // позаекранне зображення

            private Graphics m_g; // графічний об’єкт позаекранного зображення

double f(int i,  double х)  // метод класу, що вказує на функцію, графік якої          // потрібно відобразити

            {

                        return SampleSystemSolution.ylag(i, õ);

             }

Об’єкт  класу DelaySystemSolution, асоційований із системою із запізненням, що досліджується, створюється в класі GraphConstruction і є його змінною.

            Далі зупинимося на будові класу DelaySystemSolution. Нижче вказані його головні змінні.

public class DelaySystemSolution

{

            double[] у; // змінна-масив для зберігання початкового значення і шуканого

// розв”язку в точці xend

            double x0; // початкова точка для х

            private double xstor[]; // масив точок, в яких проводилося інтегрування

            private double ystor[][]; // масив знайдених розв”язків

            private double c1[][]; // масиви коефіцієнтів апрксимуючих  многочленів,

            private double c2[][]; // в яких зберігаються коефіцієнти при доданках

            private double c3[][]; // відповідно першої, другої, третьої, четвертої

            private double c4[][]; //  ñòåïåí³

            int n; // розмірність системи; може бути <= 51

Метод fcn, оголошений як

 public double[] fcn(double õ,  double y[])

повертає масив значень правих частин системи (1). Нижче буде показана його реалізація на прикладі розв”язування задачі медичної інформатики.

            Метод phi, оголошений як

            public double phi(int i,  double õ)

повертає значення початкової функції  , обчислене для i-ї компоненти розв”язку.

Знаходження наближеного розв”язку задачі (1), (2) здійснюється в класі DelaySystemSolution таким чином. У конструктор класу DelaySystemSolution включений виклик методу retard, оголошеного як:

public void retard(int n,  double õ,  double y[], double xend,  double eps,  double hmax,  double h)

Далі, після успішного кроку інтегрування метод retard викликає метод store, який записує у змінні xstor,  c1, c2, c3, c4 коефіцієнти многочленів з відповідними значеннями  .

Після цього при виклику методу  ylag(i, x) здійснюється пошук потрібного інтервалу  для   і обчислюється відповідний многочлен  для i-го розв”язку. Отже, за допомогою методу  ylag можна отримати  i-й розв”язок в довільній точці   . Ця властивість використовується для задання запізнення в методі fcn.

Як приклад розглянемо задачу про реакцію клітки на рентгенівське опромінення. Припустимо, що опромінення триває протягом часу  (і починається при  ), а потім припиняється.

Припускаємо, що клітини мають здатність поповнювати недостаток або усувати надлишок певної  речовини, але їх реакція відбувається із запізненням, що дорівнює   . Нехай   – концентрація речовини в клітині, що піддається опроміненню,  – стала опромінення, що залежить від степені опромінення,  – стала, що показує реакцію клітини на відхилення від рівноважної концентрації  . Очевидно, що  . Виходячи із зроблених припущень, концентрація описується рівняннями:

                                                            (3)

Розв”язки рівнянь можуть бути знайдені методом, запропонованим в даній роботі. Для цього потрібно створити такий метод fcn:

        public double[] fcn(double x, double y[])

        {

             //this method describes “right side” of system

             double T = 20.;

             double tau = 1.;

             double[] dRight_side = new double[Math.max(n, nn) + 1];

             if((x>0) & (x<T))

             {

                  dRight_side[1] = m_a * y[1] + m_b * (ylag(1, x – tau) – phi(1,0));

             }

             else

             {

                  dRight_side[1] = m_b * (ylag(1, x – tau) – phi(1,0));

             }

             return dRight_side;

       }

Припускаємо. що что m_a і m_b є змінними класу DelaySystemSolution, призначеними для зберігання значень  . На малюнку 2 представлено розв”язок вказаної задачі для випадку   .

            В якості наступного прикладу розглядається модель запального процесу інфекційної природи.      У загальному вигляді математична модель імунітету описана в роботі [2]. Вона універсальна і справедлива не тільки для запального процесу, але і для інфекційного зараження організму. У моделі враховуються наступні визначальні для перебігу процесу чинники:

I.                    популяція антигенів  , що розмножуються в організмі;

II.                 популяція антитілотвірних клітин  (плазмоклітин)  ;

III.               кількість антитіл (імуноглобулінів)  в організмі;

IV.              ступінь пошкодження  органа  .

Рівняння, що визначають динаміку перерахованих чинників, мають вигляд

  ,

  ,

  ,

з початковими умовами при  ;

  .

 Тут   – коефіциєнт розмноження антигена;   – коефіцієнт, що визначає ймовірність нейтралізації антигена антитілом;   – коефіцієнт, що зумовлює ймовірність зустрічі антиген-антитіло;   – коефіцієнт, обернений до часу життя плазмоклітин;   – швидкість виробництва антитіл однією плазмоклітиною;   – коефіцієнт, оберненопропорційний до часу розпаду антитіл;   – число антитіл, що вимагаються на нейтралізацію одного антигена;   – коефіцієнт, що визначає швидкість загибелі клітин за рахунок пошкоджуючої дії антигена;   – коефіцієнт, що враховує швидкість відновлення пошкодженого органу;   – фаза запізнення (час, за який здійснюється формування аскаду плазмоклітин);   – неперервна незростаюча функція ( ), що характеризує порушення нормального функціонування імунної системи через значне пошкодження органа-мішені.

Перераховані параметри додатні і є специфічними як для виду антигенна, так і для органу і конкретного організму.        За допомогою розробленої комп’ютерної моделі здійснене якісне дослідження запального процесу у випадку, коли:

  .

При   справедливі такі початкові умови

  . 

Здійснене моделювання показує, що час повторної появи запального процесу і ступінь його активності залежать від коефіцієнту  , що узгоджується з експериментальними даними (див. малюнки 3, 4).

1.     Створіть новий каталог для збереження усіх файлів проекту.

2.     Розмістіть у каталозі файли Java-класів: BoundsLocation, ColorList, DelaySystemSolution, GraphConstruction, LinkedList, ListOfSolutions, xyPoint.

3.     Створіть тектовий файл вигляду:

<html>

<head>

<title>Система імунного захисту</title>

</head>

<body>

<hr>

<applet

    code=GraphConstruction.class

    id=GraphConstruction

    width=700

    height=900 >

<param name=x0 value=0>

<param name=x1 value=60>

<param name=Scale value=30>

<param name=Legend value=true>

<param name=GraphWidth value=300>

<param name=GraphHeight value=300>

<param name=GraphCount value=4>

<param name=PunctureLine value=true>

<param name=Delay value=.5>

<param name=Hmax value=.2>

<param name=Beta value=2.>

<param name=Gamma value=.8>

<param name=Alpha value=10000>

<param name=Mu_c value=.5>

<param name=Ro value=.17>

<param name=Mu_f value=.17>

<param name=Eta value=10.>

<param name=Sigma value=300.>

<param name=Mu_m value=.12>

<param name=Tau value=.5>

</applet>

</body>

</html>

4.     Збережіть файл, надавши йому розширення .html.

5.     Запустіть файл через програму Web-броузер, встановлену в системі. Результат виконання аплету показано на малюнку

t t t t t t t t t tJava is disabled in your Browser! t

1.     Дослідіть зміну різних характеристик імунної системи у різних часових межах:

2.     Здійсніть одночасне виведення кількох графіків у одній площині за допомогою кнопки Mark few graphs:

3.     Здійсніть повноекранний показ графіка за допомогою кнопки Whole screen:

4.     Змоделюйте графіки для імунної системи для наступних значень параметрів моделі:

Примітка. Для цього вносяться відповідні зміни у речення типу

<param name=Beta value=2.>

файлу .html.

5.     Створіть Web-сторінку, яка б одночасно моделювала імунні системи для двох груп параметрів, поданих у завданні 4.

1.   Минцер О.А., Молотков В.Н., Угаров Б.Н. и др. Биологическая и медицинская кибернетика. Справочник. – К.: Наукова думка, 1989. – 375 с.

2.   Марчук Г.И. Математические модели в иммунологии. – М.: Наука, 1980. – 264с.

3.   Хайрер Д., Нерсетт В., Ваннер Т. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. – М.: Мир, 1990. – 412с.

4.   Хейл Дж. Теория функционально-дифференциальных уравнений. – М.: Мир, 1984. – 421с.

5.   Красовский Н.Н. Некоторые задачи теории устойчивости движения. – М.: Физматгиз, 1959. – 212с.

6.   Дэвис Стефен Р. Программирование на Microsoft Visual Java++/Пер.с англ. – М.: Издательский отдел “Русская редакция” ТОО “Channel Trading Ltd.”,  1997. – 376с.: ил.

Малюнок 1. Схема класів проекта розв”язування задач (1), (2).

Малюнок 2. Зміна концентрації речовини в результаті рентгенівського опромінення.

Малюнок 3. Модель запального процесу інфекційної природи 

при  . Тут x1 , x2 ,  x3 ,  x4 .

Малюнок 4. Модель запального процесу інфекційної природи 

при  . Тут x1 , x2 ,  x3 ,  x4 .