ANALISIS NUMERIK SISTEM DINAMIK DAN SINKRONISASI PROPAGASI TIPE 1 DAN 2 MODEL SARAF TERKOPEL MORRIS-LECAR ADAM SUKMA PUTRA

ANALISIS NUMERIK SISTEM DINAMIK DAN SINKRONISASI PROPAGASI TIPE 1 DAN 2 MODEL SARAF TERKOPEL MORRIS-LECAR

ADAM SUKMA PUTRA

DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011

Abstrak Adam Sukma Putra. Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar. Dibimbing oleh: Dr. Agus Kartono (Pembimbing I), Ardian Arief M.Si. (Pembimbing II).

Analisis model saraf tipe 1 dan 2 pada tugas akhir ini meliputi analisis sistem dinamik dan sinkronisasi pada variasi arus terapan DC konstan, DC dan AC bergantung waktu. Arus DC konstan menghasilkan eksitasi periodik dengan karakteristik dinamik yang tidak berubah. Pada arus bergantung waktu mengalami perubahan pola propagasi dan karakteristik dinamiknya. Pada tipe 1 merupakan bifurkasi saddle-node, pada tipe 2 andronov-hopf. Tipe 1 memiliki pita arus eksitasi lebih sempit dari tipe 2. Analisis sinkronisasi untuk sistem banyak saraf (n=2,3,4), untuk tipe 1 dan 2 meliputi kekuatan kopling εij keterhubungan saraf hij dan fase propagasi. Sinkronisasi pada kedua tipe dapat terjadi saat keadaan phase-locking tercapai saat tiap sel saraf pada sistem kopling memiliki kekuatan ikatan (ε) yang sama (εi=εj). keadaan phase-locking dibagi dua yaitu different phase synchronization jika saraf terkopel memiliki fase berbeda dan same phase synchronization jika memiliki fase yang sama.

Kata kunci: saraf, arus listrik, sistem dinamik, bifurkasi, saraf terkopel , sinkronisasi.

Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar

Adam Sukma Putra

Skripsi

Sebagai salah satu syarat untuk memperoleh gelar Sarjana Sains pada Departemen Fisika

DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011

PERNYATAAN Dengan ini penulis menyatakan bahwa skripsi yang berjudul Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar adalah karya ilmiah yang sebenar-benarnya dibuat oleh penulis sendiri dibawah bimbingan Dr. Agus Kartono dan Ardian Arief, M.si.. untuk menghindari plagiatisme, segala bentuk kutipan baik langsung maupun tidak langsung telah penulis cantumkan nama sumber maupun peneliti yang berhubungan dengan kutipan tersebut dalam bentuk teks dan dalam daftar pustaka. Skripsi ini belum pernah dipublikasikan kepada individu ataupun instansi manapun baik secara legal maupun ilegal. Bogor, April 2011

Adam Sukma Putra

Judul : Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar Nama : Adam Sukma Putra NRP

: G74070016

Disetujui,

Dr. Agus Kartono

Ardian Arief, M.Si

Pembimbing I

Pembimbing II

Diketahui,

Dr. Ir. Irzaman, M.Si. Ketua Departemen Fisika

Tanggal Lulus:

Judul : Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar Nama : Adam Sukma Putra NRP

: G74070016

Disetujui,

Dr. Agus Kartono

Ardian Arief, M.Si

Pembimbing I

Pembimbing II

Diketahui,

Dr. Husin Alatas. Kepala Bagian Fisika Teori

Tanggal Disetujui:

RIWAYAT HIDUP Penulis dilahirkan di Jakarta, pada hari Jum’at tanggal 08 Desember 1989. Penulis merupakan anak pertama dari pasangan Sukmana Nata Permana BE. dan Widayati. Penulis memiliki 2 saudara laki-laki yaitu Gilang Sukma Putra (Mahasiswa) dan Adytia Gumelar (Pelajar). Riwayat pendidikan penulis di mulai pada tahun 1995 di SDN Cibatok I Bogor selama enam tahun dan lulus berijazah pada tahun 2001. Kemudian penulis meneruskan ke SMPN I Cibungbulang selama tiga tahun dan lulus tahun 2004 kemudian masuk SMAN I Leuwiliang dan lulus pada tahun 2007. Penulis kemudian melanjutkan studi di Departemen Fisika FMIPA IPB hingga tahun 2011. Selama penulis kuliah di IPB, penulis aktif mengikuti organisasi kemahasiswaan diantaranya Himpunan mahasiswa Fisiska (HIMAFI) sebagai staf instrumentasi dan Teknologi (INSTEK) tahun kepengurusan 2009, dan hingga saat ini penulis aktif di UKM Uni Konservasi Fauna (UKF) IPB dan pernah menjabat sebagai ketua divisi PSDM, ketua Espedisi Global UKF 2008, dan ketua kegiatan keilmuan di UKF. Penulis juga aktif sebagai pemateri dan moderator pada berbagai kegiatan keilmuan antara lain sebagai moderator pelatihan PKM. Penulis juga produktif dalam prestasi akademik antara lain pernah menyandang Mahasiswa Berprestasi (MAPRES) Departemen Fisika periode 2009-2010, dan aktif sebagai peserta OSN-PTI tahap I dan II. Saat ini, penulis terus berusaha untuk mengembangkan ilmu fisika dalam segala terapan kehidupan terutama pada Bagian Fisika Komputasi dan Pemodelan.

DAFTAR ISI Halaman DAFTAR GAMBAR .................................................................................................... x DAFTAR TABEL ...................................................................................................... xiii DAFTAR LAMPIRAN .............................................................................................. xiv I

PENDAHULLUAN......................................................................................... 1 1.1 Latar Belakang........................................................................................... 1 1.2 Tujuan Penelitian ....................................................................................... 2 1.3 Perumusan Masalah ................................................................................... 2 1.4 Hipotesis Penelitian ................................................................................... 2

2

TINJAUAN PUSTAKA .................................................................................. 2 2.1 Saraf ........................................................................................................... 2 2.1.1 Neurofisiologi ................................................................................... 2 2.1.2 Mekanisme Ionic Transport pada saraf ............................................ 3 2.2 Model Morris-Lecar (1981) ....................................................................... 5 2.3 Sistem Dinamik dan Bifurkasi ................................................................... 6 2.3.1 Equilibrium ....................................................................................... 6 2.3.2 Analisis linier lokal ........................................................................... 7 2.3.3 Nilai eigen dan vektor eigen ............................................................. 7 2.3.4 Klasifikasi equilibrium ..................................................................... 8 2.4 Propagasi Saraf .......................................................................................... 9 2.4.1 Bifurkasi ......................................................................................... 10 2.4.2 Klasifikasi propagasi saraf.............................................................. 11 2.5 Sinkronisasi ............................................................................................. 12

3

METODE PENELITIAN .............................................................................. 14 3.1 Waktu dan Tempat Penelitian.................................................................. 14 3.2 Peralatan .................................................................................................. 14 3.3 Metode Penelitian .................................................................................... 14

ix

3.3.1 Studi pustaka................................................................................... 14 3.3.2 Analisis numerik solusi propagasi saraf Morris-Lecar dengan metode RK-4 ................................................................................. 14 3.3.3 Analisis sistem dinamik model saraf Morris-Lecar ........................ 14 3.3.2 Sinkronisasi model saraf Morris-Lecar terkopel .......................... 14 4

HASIL DAN PEMBAHASAN ..................................................................... 15 4.1 Solusi Numerk Propagasi Saraf dengan Metode RK-4 ........................... 15 4.1.1 Solusi numerik dengan arus terapan DC tetap................................ 16 4.1.2 Solusi numerik dengan arus terapan DC bergantung waktu........... 17 4.1.3 Solusi numerik dengan arus terapan AC bergantung waktu........... 20 4.2 Analisis Sistem Dinamik Propagasi Saraf ............................................... 21 4.2.1 Analsis linier lokal, nilai eigen, dan diagram fase .......................... 22 4.2.2 Nilai eigen dan diagram fase tipe 1 dan 2 variasi Iapp dan Vc ......... 25 4.2.3 Nilai eigen dan diagram fase tipe 1 dan 2 Iapp bergantung waktu ........................................................................ 26 4.3 Solusi Numerik Propagasi Saraf Terkopel .............................................. 30 4.3.1 Model saraf terkopel ....................................................................... 31 4.3.2 Solusi numerik model saraf terkopel Iapp tetap ............................... 32 4.3.3 Solusi numerik model saraf terkopel Iapp AC bergantung waktu ... 35 4.4 Solusi Numerik pada n Saraf Terkopel.................................................... 37 4.4.1 Solusi numerik pada 3 saraf terkopel ............................................. 37 4.4.2 Solusi numerik pada 4 saraf terkopel ............................................. 40

5

SIMPULAN DAN SARAN........................................................................... 43 5.1 Simpulan .................................................................................................. 43 5.2 Saran ........................................................................................................ 43

DAFTAR PUSTAKA ................................................................................................. 44 LAMPIRAN ................................................................................................................ 46 NOMENKLATUR ...................................................................................................... 61

DAFTAR TABEL Halaman 1 Hubungan nilai V3 dan Iapp dengan bifurkasi .......................................................... 26 2 Nilai eigen masing-masing daerah pada tipe 1 dan 2 arus DC bergantung waktu . 27 3 Nilai eigen masing-masing daerah pada tipe 1 dan 2 arus AC bergantung waktu . 29

DAFTAR GAMBAR Halaman 1 Struktur saraf ........................................................................................................... 3 2 Sinapsis ................................................................................................................... 3 3 Nilai Nernst Potential tiap ion di lingkungan dalam dan luar membran pada mamalia saat T=370C (310 K) ............................................................................... 4 4 Mekanisme terjadinya potensial aksi pada saraf ..................................................... 5 5 Titik kritis pada keadaan Stabil ............................................................................... 7 6 Titik kritis pada keadaan tak Stabil. (a) tiga equilibrium -66 mV, -56 mV, dan -28 mV (b) satu equilibrium -61 mV ............................................................................. 7 7 Klasifikasi titik equilibrium .................................................................................... 8 8 Tititk equilibrium node ,λ1=-1, λ2=-3 (stabil), λ1=+1, λ2=+3 (tidak stabil) ............ 8 9

Saddle equilibrium λ1=+1, λ2=-1 ............................................................................ 9

10 Equilibrium λ1=-3±i (stabil) atau λ1=+3±i (tidak stabil) ......................................... 9 11 Suatu limit cycle. (a) keadaan terangsang Excitability. Sebuah trayektori (kotak kecil) Terinisiasi dekat dengan titik kestabilan dan kembali ke tiik kestabilan tersebut. (b) sistem tereksitasi dekat dengan bifurkasi akan membentuk limit cycle karena titik keseimbangan nya belum terinisiasi (selama stimulus diterapkan secara kontinu) ...................................................................................................... 10 12 Transisi terjadi dari keadaan istirahat menuju keadaan eksitasi yang terjadi melalui bifurkasi keseimbangan (tanda panah). Bifurkasi Saddle-Node pada in vitro pyramidal neuron of rat’s primary visual cortex. Bifurkasi Andronov-Hopf pada in vitro brainstem mesencephalic V neuron. ................................................ 10 13 Limit cycle empat bifurkasi pada keadaan istirahat............................................... 11 14 Sistem dinamik tereksitasi dapat memiliki bifurkasi yang berbeda baik melalui keadaan bistable atau tidak ................................................................................... 11 15 Propagasi eksitasi saraf (a) tipe 1 dan (b) tipe 2 dengan arus terapan DC konstan pada sela saraf tikus .............................................................................................. 12 16 Tiga kondisi sinkronisasi (a) sefase (b) berlawanan fase (c) berbeda fase ........... 13

xi

17 Dua saraf terkopel dengan peta frekuensi impuls (a) berbeda fase (b) sefase (c) berlawanan fase ..................................................................................................... 17 18 Aktivitas listrik (action potential) model saraf Morris-Lecar tipe 1 ..................... 16 19 Bentuk propagasi saraf tipe 2 ................................................................................ 17 20 Nilai Iapp pada (a) tipe 1 dan (b) tipe 2 masing-masing 40 mV dan 50 mV. Kedua bentuk propagasi tidak dapat terjadi secara periodik ............................................ 17 21 Propagasi saraf tipe 1 dengan arus I(t). ................................................................. 18 22 Propagasi saraf tipe 2 dengan arus I(t). ................................................................. 18 23 Frekuensi propagasi (spike/second) pada (a) tipe 1 dan (b) tipe 2 ........................ 19 24 Propagasi (a) tipe 1 dan (b) tipe 2 dengan gradient I(t) negatif ........................... 19 25 Propagasi saraf dengan fungsi arus terapan AC.(a) tipe 1.(b) tipe 2..................... 20 26 Variasi nilai ω terhadap bentuk propagasi saraf ................................................... 21 27 Diagram fase (a) tipe 1 dan (b) tipe 2 dengan Iapp tetap ........................................ 22 28 Bifurkasi saddle-node pada tipe 1 ......................................................................... 24 29 Bifurkasi saddle-node on invariant circle (SNIC) ................................................ 24 30 Bifurkasi Andronov-Hopf pada tipe 2 .................................................................. 24 31 Bifurkasi Andronov-Hopf ...................................................................................... 25 32 Diagram bifurkasi (a) tipe 1 dan (b) tipe 2 dengan nilai V3=18 mV ..................... 25 33 Tiga daerah utama propagasi (a) tipe 1 dan (b) tipe 2 dengan arus DC bergantung waktu ..................................................................................................................... 26 34 Diagram fase tipe 1 dengan fungsi arus DC bergantung waktu ............................ 28 35 Diagram fase tipe 2 dengan fungsi arus DC bergantung waktu ............................ 28 36 Diagram fase tipe 1 dengan fungsi arus AC bergantung waktu ............................ 29 37 Diagram fase tipe 2 dengan fungsi arus AC bergantung waktu ............................ 30 38 Dua saraf tipe 1 non-kopling. εi,2=0 ...................................................................... 32 39 Tipe 1 dua saraf terkopel (εi,2 ≠ 0, hij=1), non-sinkronisasi .................................. 32 40 Tipe 1 dua saraf terkopel (εi,2 ≠ 0, hij=1), tersinkronisasi (εi = ε2=0.5 mS/cm2) ... 33 41 Sinkronisasi kopling saraf tipe 1 dengan nilai skala waktu berbeda ..................... 33 42 Tipe 2 dua saraf terkopel (εi,2 ≠ 0, hij=1), non-sinkronisasi .................................. 34 43 Tipe 2 dua saraf terkopel (εi,2 ≠ 0, hij=1), tersinkronisasi (εi = ε2=0.5 mS/cm2) .. 34

xii

44 Sinkronisasi kopling saraf tipe 2 dengan nilai skala waktu berbeda ..................... 35 45 Sinkronisasi tipe 1 dengan arus terapan AC ......................................................... 36 46 Sinkronisasi tipe 2 dengan arus terapan AC. ........................................................ 36 47 Sistem 3 saraf terkopel. (a) terkopel dan sinkron (b) kopel tidak sempurna......... 38 48 Sinkronisasi 3 saraf terkopel tipe 1 dengan variasi kemungkinan keadaan terkopel. (a) arus terapan AC. (b) tidak terkopel. (c) terkopel dengan fase berbeda. (d) terkopel dan tersinkronisasi. (e) kopel tidak sempurna......................................... 38 49 Sistem 3 saraf terkopel. (a) variasi 1. (b) variasi 2................................................ 39 50 Sinkronisasi 3 saraf terkopel tipe 1 dengan variasi kemungkinan keadaan terkopel. (a) variasi 1 beda fase (b) variasi 1 sefase (c) varaisi 2 beda fase (d) varaisi 2 sefase ..................................................................................................................... 39 51 Sinkronisasi 3 saraf terkopel tipe 2 dengan variasi kemungkinan keadaan terkopel. (a) arus terapan AC (b) tidak terkopel (c) terkopel dengan fase berbeda (d) terkopel dan tersinkronisasi (e) kopel tidak sempurna .......................................... 40 52 Sinkronisasi 3 saraf terkopel tipe 2 dengan variasi kemungkinan keadaan terkopel (a) variasi 2 beda fase (b) variasi 2 sefase............................................................ 40 53 Sistem 4 saraf terkopel .......................................................................................... 41 54 Sinkronisasi 4 saraf terkopel tipe 1 dan 2 (a) tipe 1 beda fase (b) tipe 1 sefase (c) tipe 2 beda fase (d) tipe 2 sefase ........................................................................... 42 55 Variasi sistem 4 saraf terkopel .............................................................................. 42 56 Variasi propagasi 4 saraf terkopel (a) tipe 1 (b) tipe 2 .......................................... 42

DAFTAR LAMPIRAN Halaman 1 Diagram alir Penelitian ......................................................................................... 47 2 Metode Rungge-Kutta Orde 4 (RK-4) .................................................................. 48 3 Metode Rungge-Kutta orde 4 untuk model Morris-Lecar (1981) ......................... 49 4 Program untuk mencari solusi numerik propagasi saraf ....................................... 51 5 Program untuk mencari nilai akar keseimbangan ................................................. 52 6 Program untuk mencari matriks jacobian dan nilai eigen ..................................... 53 7 Program untuk grafik garis nol V nulclines dan W nulclines ................................ 55 8 Program untuk saraf terkopel n=2,3,4 ................................................................... 55 9 Program untuk saraf 1,2,3, dan 4 .......................................................................... 58

KATA PENGANTAR Dengan segala kekuasaan Allah S.W.T dan maghfirah Nabi Muhammad SAW., akhirnya penulis dapat menyelesaikan skripsi dengan judul “Analisis Numerik Sistem Dinamik dan Sinkronisasi Propagasi Tipe 1 dan 2 Model Saraf Terkopel Morris-Lecar” oleh karena itu sebesar-besar nya penulis panjatkan puji syukur ke hadirat Nya atas nikmat yang besar ini. Skripsi ini memiliki topik yang menarik untuk dikaji karena berhubungan dengan sistem informasi saraf manusia yang dimodelkan secara matematis dan fisis. sehingga merupakan suatu ilmu yang sangat perlu dikembangkan terus dalam kanca perkembangan sains di Indonesia. Penulis berharap, penelitian ini dapat bermanfaat bagi perkembangan pengetahuan di Indonesia terutama pada bidang fisika teori dan terapannya. Peneliti ucapkan terima kasih yang sebesar-besarnya atas semua pihak yang terlibat dalam kelancaran tugas akhir ini, yaitu kepada: 1.

Ayahanda dan Ibunda tercinta yang tidak pernah habis kasih sayang menjaga, merawat, dan memberi pendampingan moril maupun materil selama penulis hidup.

2.

Dr. Agus Kartono selaku pembimbing utama yang dengan sabar dan bijaksana dalam membimbing penulis dalam menyelesaikan tugas akhir ini.

3.

Ardian Arief Setiawan, M.Si. selaku pembimbing kedua dengan semangat dan dorongan yang telah diberikan kepada penulis.

4.

Hanedi Darmasetiawan dengan kesabaran dan ketelitiannya membantu penulis dalam merampungkan penyusunan skripsi ini.

5.

Saudara tercinta Gilang Sukma Putra dan adinda tercinta Adytia Gumelar atas senyuman kasih sayang dan semangat yang tak pernah terputus selama penulis hidup.

6.

Teman-teman satu bidang fisika teori Departemen Fisika IPB Subiyanto, Putra Wira Kurniawan, Dede Hermanudin, Ika Ikrima, Nurullaeli, Switenia Wana Putri, dan Izzatu Yazidah dan seluruh teman seangkatan 44 tersayang atas

vii

dukungan spirituil dan hangat persahabatan yang telah peneliti rasakan. Adik angkatan 45 dan 46 seperjuangan, semoga cepat lulus. 7.

Seluruh anggota Unit Kegiatan Mahasiswa Uni Konservasi Fauna (UKF) IPB yang telah memberikan kenyamanan dalam menerapkan ilmu konservasi di dalam maupun di luar kampus.

8.

Teman-teman penerima beasiswa Tanoto Foundation yang telah memberikan semangat moril maupun materil selama penulis menyelesaikan studi di IPB.

9.

Seluruh staf Dosen pengajar Departemen Fisika IPB yang telah membuat penulis kagum dan menikmati indahnya fisika.

10. Seluruh staf karyawan Departemen Fisika IPB yang telah membuatan penulis nyaman di fisika. 11. Seluruh staf yang telah sangat membantu kelancaran tugas akhir ini yang tidak dapat penulis sebutkan satu-persatu. Ketidaksempurnaan selalu melekat pada skripsi yang telah penulis buat ini, oleh karena itu segala masukan dan kritik sangat diharapkan bagi penulis guna menjadikan skripsi ini lebih baik. Akhir kata, semoga kelak nantinya penelitian ini dapat bermanfaat dan menjadi terobosan baru bagi perkembangan sains di Indonesia.

Bogor, 28 April 2010

penulis

1

I. PENDAHULUAN 1.1 Latar Belakang Sistem saraf merupakan sistem utama yang berperan dalam pengaturan mekanisme kerja tubuh. Mekanisme kerja tubuh tersebut terkoordinasi dalam suatu sistem pusat yang membentuk fungsifungsi spesifik, yang disebut fungsi tubuh.1 Berbagai macam stimulus baik dari luar maupun dari dalam tubuh terkoordinasi dan terintegrasi sehingga membentuk respon aktif maupun pasif sebagai umpan balik tubuh. Stimulus ini berupa impuls-impuls saraf yang menjalar akibat adanya perbedaan nilai membran potensial pada sel saraf.1 Berbagai fenomena ini sangat kompleks dalam sistem saraf manusia dan merupakan suatu hal yang menarik untuk dipelajari. Neuron atau dikenal dengan nama sel saraf meupakan satuan unit dari sistem saraf yang membangun tubuh manusia.1 Saraf memiliki peranan penting dalam sistem koordinasi. Saraf merupakan dasar dari pembangunan sistem saraf pusat yang bertanggung jawab atas kontrol gerak tubuh, sistem kesadaran, ingatan, hingga intelektualitas manusia.1 Meskipun pengkajian yang telah dilakukan untuk memahami bagaimana sistem saraf berinteraksi dalam pembentukan daya inelektualitas sangat terbatas, namun salah satu cara yang sederhana untuk memahaminya adalah dengan cara memodelkan jaringan saraf dengan mempelajari karakteristik penjalaran pada sel saraf individu.2 Jaringan saraf telah banyak dipelajari oleh para peneliti karena sangat menarik dan mudah dalam penerapannya. Sebagaimana yang telah disebutkan sebelumnya, bahwa penjalaran impuls saraf melibatkan perubahan potensial membran yang disebabkan oleh adanya aktivitas transpor ionik pada synapses (tempat terjadinya reaksi elektrokimia saat impuls menjalar dari satu sel ke sel lain) pada tiap sel nya. Dalam sistem dinamik saraf, para peneliti memodelkan mekanisme penjalaran impuls saraf sebagai suatu gangguan periodik yang

spontan (spontaneous spiking behavior) karena adanya perubahan membran potensial selama penjalaran.2,3 Fenomena ini dapat digambarkan sebagai kondisi ketika tidak ada stimulus maka tidak akan ada penjalaran impuls atau saraf dalam keadaan istirahat (resting point). Kejadian pada periode ini sama pada setiap sel saraf. Oleh karena itu, untuk mempelajari penjalaran pada jaringan saraf, dapat direpresentasikan oleh satu model sel saraf.1,2 Mekanisme penjalaran impuls pada saraf sangat erat kaitannya dengan perubahan beda potensial pada saraf, energi kinetik dari arus, potensial aktivasi dan istirahat, serta parameter lainnya. Para peneliti telah banyak membuat model saraf ini dengan meninjau mekanisme dinamik yang terjadi pada saraf. Seperti Boltzman activation function yang meninjau perubahan parameter arus pada membran saraf , fungsi laju perubahan potensial membran , dan nilai konduktansi membran.1 Model saraf ini dibangun dalam sebuah persamaan matematis yang dapat dianalisis secara numerik.3 Pemodelan pertama saraf dalam komputasi neuroscience adalah model Hodgkin-Huxley (HH) (1952) yang biasa disebut The Squid Giant Axon Model (SGA). Model ini dibangun dari hasil eksperimen pada gurita raksasa. Hodgkin dan Huxley menyatakan bahwa dalam penjalarannya, SGA membawa tiga arus utama pada mebran saraf yaitu: parameter yang berkaitan dengan aktivasi K+, parameter inaktivasi Na+, dan parameter arus reset yang dikendalikan oleh Cl-. Persamaan ini terangkum dalam sistem empat dimensi dari perubahan membran potensial, dan laju ionik tiap parameter aktivasi. Selain model ini, model lain yang cukup pupuler adalah model Fitzhugh-Nagumo (FN) yang merupakan sebuah model penyederhanaan dari HH. Model FN disajikan dalam dua variabel utama yaitu variabel pemercepat (excitation variable) dan variabel pelambat (recovery 3 variable).

2

Pada penelitian ini, digunakan model yang telah di pubikasikan oleh Cathy Morris dan Harold Lecar (1981) (ML Model).24 Model ini diturunkan dari hasil eksperimen mengenai sifat listrik dari serat otot angsa putih (Barnacle Muscle Fiber (BMA)) yang menunjukan aktivitas listrik saat diterapkan suatu arus luar. Serat ini terutama mengandung beda potensial saluran kalium (voltage gated K+) dan arus Ca+ dengan sebuah arus K+ yang diaktifkan oleh Ca+ di bagian dalam sel. Kedua parameter ini sangat berperan dalam aktivasi potensial listrik pada BMA yang dapat terjadi melalui suatu mekanisme perbedaan konduktansi Ca+ dalam aktivasi dan beda potensial K+ untuk pemulihan.24 Model ini terdiri dari persamaan dua dimensi yang melibatkan parameter arus dan kapasitansi membran, dengan parameter aktivasi utama berdasarkan nilai konduktansi dari saluran kalsium dan kalium pada membran. Secara matematis, pada model ini digunakan suatu fungsi persamaan kemungkinan yang diturunkan dengan asumsi bahwa dalam sebuah kesetimbangan (keadaan stabil), keadaan terbuka dan tertutup nya saluran pada membran dibatasi distribusi Boltzmann. Dengan demikian, keterlibatan parameter dalam model dua dimensi ini dapat memudahkan dalam analisisnya baik secara fisika maupun mengenai sistem dinamiknya terhadap interpretasi fungsi biologi dari membran pada saraf dalam suatu jaringan tubuh. Metode analisis yang dipakai pada penelitian ini adalah dengan penerapan penggunaan perangkat lunak untuk komputasi. Metode yang digunakan adalah analisis numerik persamaan differensial (PD) yang terbagi atas visualisasi penjalaran impuls dengan menggunakan metode Rungge-Kutta 4 (RK-4), analisis sistem dinamik meliputi bifurkasi dan ruang fase, serta metode sinkronisasi penjalaran impuls pada saraf kompleks (lebih dari satu sel saraf). Simulasi yang dilakukan pada penelitian ini menggunakan MATLAB sebagai sarana pengolah data dari analisis kuantitatif model persamaan matematis.

1.2 Tujuan Penelitian Penelitian ini bertujuan untuk menganalisis pengaruh arus listrik terapan tetap serta arus DC dan AC bergantung waktu pada propagasi saraf tipe 1 dan 2. Analisis sistem dinamik saraf tipe 1 dan 2 pada model MorrisLecar (1981), serta sinkronisasi sistem saraf terkopel pada jaringan kompleks (lebih dari satu saraf). 1.3 Perumusan Masalah 1. Bagaimana tingkat ketepatan metode RK-4 terhadap pola propagasi saraf tipe 1 dan 2 dengan berbagai variasi nilai arus terapan pada model? 2. Apakah variasi nilai arus terapan pada model Morris-Lecar mempengaruhi pola potensial aksi pada saraf ? 3. Bagaimana karakteristik sistem dnamik pada model Morris-Lecar dapat menjelaskan secara kualitatif pola propagasi saraf tipe 1 dan 2? 4. Bagaimana pola propagasi dan sinkronisasi saraf tipe 1 dan 2 pada sistem saraf kompleks terkopel? 1.4 Hipotesis 1. Metode RK-4 memiliki keakuratan yang tinggi pada simulasi pola penjalaran saraf. 2. Karakteristik sistem dinamik pada model ML untuk tipe 1 adalah saddle-node dan untuk tipe 2 adalah Andronov-Hopf. 3. Pola propagasi pada sistem kompleks memiliki jenis sinkronisasi sefase dan berbeda fase pada tiap sel saraf terkopel.

II. TINJAUAN PUSTAKA 2.1 Saraf 2.1.1 Neurofisiologi Saraf “neuron” adalah suatu sel saraf yang merupakan unit anatomis dan fungsional dari sistem saraf. Secara fisiologi,saraf terdiri dari tiga bagaian utama yaitu: dendrit, badan sel atau soma, dan axon.1 Dendrit merupakan bagian masukan dari sebuah saraf dan penerima masukan sinaptik dari luar

3

lingkungan saraf atau dari saraf lain.2 Bagian soma merupakan konstruksi utama pada saraf yang didalamnya terdapat unit-unit seluler yang penting dan seperti inti sel (nucleus) mitokondria.2 Bentuk perpanjangan badan sel yang panjang dan berseludang disebut akson yang berfungsi sebagai penghantar informasi ke luar sel menuju lingkungan luar sel atau sel lain.1 Terkadang ada beberapa sinyal dapat diterima oleh akson.2 Struktur saraf secara umum dapat dilihat pada Gambar 1.

efektor.1 Melalui mekanisme difusi listrik, pada synapses ini saraf satu dengan yang lain saling terhubung.

Gambar 2. Sinapsis.1

Gambar 1. Struktur Saraf.3 Ketiga komponen utama ini saling bersinergi sehingga memiliki kemampuan spesifik untuk menerima, meneruskan, dan menyampaikan pesan dalam bentuk neural impuls karena sifat khas saraf yang sangat peka terhadap rangsangan dan dapat menghantarkan pesan elektrokimia.2 Pada akson bagian terminal terdapat synapses, yaitu bagian akson yang telah mengalami spesifikasi struktur tempat neurotransmitter dilepaskan untuk berkomunikasi dengan saraf lainnya. Synapses merupakan satusatunya tempat suatu impuls dapat lewat dari suatu saraf ke saraf lainnya atau

2.1.2 Mekanisme Ionic Transport pada saraf Aktivitas listrik yang terjadi pada saraf dalam penjalarannya melibatkan arus ionik pada membran saraf antara lain Na+, K+, Ca+, dan Cl-.3 Konsentrasi tiap ion berbeda pada lingkungan di dalam dan di luar membran. Pada bagian luar membran (ekstraseluler) memiliki konsentrasi ion Na+ dan Cl- yang lebih besar.Sedangkan pada bagian dalam membrane (intraseluler) memiliki + konsentrasi K yang tinggi dan molekul bermuatan yang bersifat negatif (A-). Pada keadaan setimbang, bagian dalam maupun luar membran memiliki konsentrasi dan nilai potensial keseimbangan tertentu saat tidak menerima rangsangan dari luar.Kondisi ini biasa disebut dengan Nernst Potential.3 Nilai-nilai potensial keseimbangan tiap ion dapat dilihat pada Gambar 3.

4

Gambar 3.Nilai Nernst Potential tiap ion di lingkungan dalam dan luar membran pada mamalia saat T=370C (310 K) .3 Mekanisme penjalaran saraf terkait dengan mekanisme transport ionik dari dan ke luar membran melalui saluran ionik (ionic channelIs). Saluran ini berupa pori yang mampu melewatkan molekul pada membran. Terdapat beberapa tipe dari saluran ionik. Namun yang paling berpengaruh pada nilai potensial membran adalah saluran Na+ dam K+.4 yang berperan secara aktif dalam mekanisme penjalaran pulsa pada saraf akibat adanya perbedaan beda potensial di dalam dan di luar membran. Saat saraf dalam keadaan tidak menerima maupun mengirimkan signal disebut dalam keadaan “at rest” (istirahat).4 Pada keadaan ini lingkungan intraseluler lebih negatif. Pada keadaan at rest, ion kalium (K+) dapat melewati membran dengan mudah sedangkan ion natrium (Na+) sulit melewati membran. Begitu pula molekul-molekut negatif (A-) tidak dapat melewati membran. Tiap-tiap ion dan molekul di lingkungan membran memiliki nilai potensial, sehingga pada keadaan istirahat (setimbang), nilai beda potensial keseluruhannya dapat dikalkulasikan yaitu sekitar -70 mV

disebut resting potential (RP).4 Ini berarti bahwa saat keadaan istirahat, lingkungan intraseluler memiliki nilai potensial lebih rendah 70 mV dibandingkan lingkungan ekstraseluler. Nilai RP Ini bervariasi sebagai suatu konsekuensi saat bahwa keadaan saraf satu dengan yang lain berbeda-beda. Sebagai contoh, dalam jurnal E.M Izhikevich (2003) menyebutkan nialai RP berkisar antara 40 mV hingga -60 mV.5 Keadaan keseimbangan dapat tercapai pada saat gradien konsentrasi dan gradien potensial listrik sebanding dan berlawanan arah satu sama lain. Nilai bersih dari arus pada membran adalah nol. Keadaan ini disebut equilibrium potential yang bergantung pada jenis ion spesifik yang di sajikan dalam bentuk Nernst equation berikut;4 =

RT [Ion] ln zF [Ion]

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (1)

[Ion]indan [Ion]out adalah konsentrasi di dalam dan luar membran, Z adalah nilai valensi ion [Z(K+)= +1, Z(Na+)= -1], R adalah konstanta gas, T merupakan temperatur absolut, dan F adalah konstanta Faraday. Kondisi ketika saraf mengirimkan sinyal atau informasi terjadi apabila ada action potential (AP). Keadaan ini biasa disebut “spike” saat terjadi AP.3 AP ini terjadi akibat adanya aktivitas depolarisasi arus pada membran saraf. Ini berarti bahwa saat RP bergerak naik menuju 0 mV. Ketika nilai potential kira-kira mencapai -55 mV, saraf akan terangsang untuk menghasilkan AP. Nilai potensial ketika saraf mulai mengirimkan informasi (spiking) disebut nilai potensial ambang (threshold potential).4 Rangsangan pertama saat terjadinya AP disebabkan oleh saluran

5

natrium yang terbuka. Karena pada saat keadaan istirahat lebih banyak ion natrium di luar membran, maka ion ini masuk ke dalam lingkungan intraseluler melalui membran. Ion natrium memiliki muatan positif, menyebabkan saraf lebih bermuatan positif sehingga mengakibatkan depolarisasi. Saat setelah terjadi depolarisasi, saluran kalium mulai terbuka tetapi dengan laju yang lebih lambat dari pembukaan saluran natrium. Saat saluran kalium terbuka, ion kalium akan keluar sel. Ini menyebabkan penghambatan pada proses depolarisasi. Saat pembukaan saluran kalium mencapai maksimum, maka saluran sodium mulai tertutup sehingga menyebabkan nilai AP kembali menuju -70 mV.25

Gambar 4. Mekanisme terjadinya potensial aksi pada saraf .1 2.2Model Morris-Lecar (1981) Cathy Morris dan Harold Lecar mengusulkan sebuah model untuk menjelaskan mekanisme sifat listrik pada serat otot angsa pada tahun 1981.Model ini merupakan persamaan dua dimensi yang hanya melibatkan sebuah arus pengaktifan Ca+, sebuah arus penyearah

K+ untuk pemulihan, dan arus pasif kebocoran pada membrane (passive leak). Persamaan model ini adalah sebagai berikut. )( − ) − !( − ) − " ( − " ) + $ %% (2) !∞ ( ) − ! ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (3) !′ = '( ( ) ′

=−

∞(

Parameter V merupakan membran potensial, dan W merupakan parameter pemulihan yang merepresentasikan nilai variasi normalisasi konduktansi ion K+. parameter W sebanding dengan nilai pengaktifan kemungkinan bahwa saluran ion K+ pada keadaan terbuka (konduksi). Persamaan (3) menggambarkan proses pemulihan oleh saluran protein saat terjadi transisi antara keadaan ion konduksi dengan ion non-konduksi. Kunci utama dari eksitasi listrk yang menyebabkan action potential adalah energi dan tingkat transisi untuk proses pembukaan saluran adalah sangat bergantung pada beda potensial membran. Secara matematis, fungsi kemungkinan pembukaan saluran M∞(V) dan W∞(V) diturunkan dengan asumsi bahwa pada keadaan setimbang, pembukaan dan penutupan sebuah saluran dibatasi berdasarkan distribusi Boltzman. Fungsi konduksi ini diberikan sebagai berikut. ∞( ) =

!∞ ( ) =

1 + tanh 2

1 + tanh -

Persamaan (4) disederhanakan

2

−

/

−

3

. 2

0

0

∙∙∙∙∙∙∙∙∙∙∙∙ (4)

∙∙∙∙∙∙∙∙∙∙∙∙ (5)

dan (5) dengan

dapat cara

6

menghilangkan fungsi hiperbolik tanh menjadi suatu fungsi eksponensial yang lebih sederhana seperti pada persamaan berikut. ∞(

)=

!∞ ( ) =

1

1 + exp 8−2 1

1 + exp 8−2 -

−

/

−

3

. 2

09

09

∙∙∙∙ (6)

∙∙∙∙∙ (7)

Dalam penelitian ini digunakan fungsi seperti pada persamaan (4) dan (5). Konstanta waktu untuk + pemulihan saluran K dalam pengaruh perubahan beda potensial bergantung pada beda potensial membran. '( ( ) =

1

− ∅ cosh - 2 2 0 3

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (8)

Parameter ø merupakan skala waktu untuk proses pemulihan. Nilai ø dapat divariasikan untuk berbagai sel yang berbeda-beda dan sangat sensitif terhadap suhu lingkungan membran.26 Model ini sangat sederhana dalam menjelaskan mekanisme listrik pada membran saraf. Model propagasi saraf yang bergantung pada tiga arus ionik: ICa, merupakan penyebab utama eksitasi listrik, IK, arus utama yang berperan dalam proses pemulihan, dan IL merupakan nilai arus kebocoran membran termasuk didalamnya nilai Resting Potential. Berbagai sistem dan fenomena eksitasi potensial dapat dimodelkan dengan memvariasikan nilai konduktansi membran (gCa, gK, dan gL).24 Karena model ini berdasarkan atas konduktansi dan arus pada membran, maka baik secara teori maupun eksperimen dapat disinkronkan untuk didemonstrasikan.

2.3 Sistem Dinamik dan Bifurkasi Konsep mengenai ruang fase,titik kritis,serta Stabilitas merupakan hal yang fundamental dalam dinamika sistem. Konsep dinamika sistem ini dapat digambarkan oleh suatu set persamaan autonomous. Ini berarti suatu set persamaan yang di dalamnya tidak memiliki hubungan ketergantungan.7 Pada persamaan saraf ini yang dimaksud autonomous berarti laju variabel terkait baik potensial membran V maupun parameter pemulihan W tidak bergantung pada skala waktu untuk nilai arus terapan tetap.3 Pada model saraf Morris-Lecar, variabel dimensional yang terkait adalah potensial membran V dan parameter pemulihan W. jika V dan W di plot pada suatu bidang dua dimensi maka disebut sebagai bidang fase atau ruang fase (phase portrait) sedangkan kurva yang terbentuk merupakan trayektori bagi PDB V dan W.7 2.3.1 Equilibrium Langkah penting dalam analisis sistem dinamik adalah menentukan nilai keseimbangannya (equilibrium) yaitu pada kondisi: @(A, C) = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (9) (A, C) = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (10)

titik (x,y) adalah sebuah equilibrium.3 Penginisiasian pada kondisi awal (x0,y0) saat x’= 0 dan y’=0, dan trayektori yang terjadi tetap selama kondisi equilibrium, maka x(t)=x0, dan y(t)= y0 untuk t ≥ 0. Sifat dari trayektori ini dapat bersifat divergen atau konvergen dari titik equilibriumnya bergantung pada kestabilannya.3 Sebagai contoh, pada model HH terkait pada saluran ion K+ pada kurva I-V memiliki tiga titik nol.

7

memiliki titik equilibrium (x0,y0). Fungsi nonlinier f dan g dapat dilinierisasi dekat equilibrium sebagai berikut.

Gambar 5.Titik kritis pada keadaan stabil.3

@(A, C) = G(A − AH ) + I(C − CH ) + ℎK ℎLMNOM ∙∙∙∙∙∙∙ (13) (A, C) = P(A − AH ) + N(C − CH ) + ℎK ℎLMNOM ∙∙∙∙∙∙∙ (14)

high order dapat berupa (x-x0)2, (x-x0)(y-y0), (x-x0)3, dan seterusnya. a, b, c, dan d adalah suatu operator sebagai berkut: Q@ (A , C ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (15. G) QA H H Q@ (A , C ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (15. I) I= QC H H Q (A , C ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (15. P) P= QA H H Q (A , C ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (15. N) N= QC H H

G= Gambar 6.Titik kritis pada keadaan tak Stabil. (a) tiga equilibrium -66 mV, -56 mV, dan -28 mV (b) satu equilibrium -61 mV.3 Kedaan trayektori pada equilibrium berkaitan dengan stabilitas sistem dinamik. Sebuah equilibrium dikatakan stabil apabila setiap trayektori mendekati titik equilibrium pada t ≥ 0. Ini berarti trayektori bersifat convergen terhadap equilibrium (Gambar 5) untuk t→∞. Sebaliknya sebuah equilibrium dikatakan tidak stabil apabila trayektori bersifat divergen atau menyebar dari equilibrium (Gambar 6.) 2.3.2 Analisis linier lokal Agar lebih memahami mengenai analisis sistem dinamik, perlu diketahui karakteristik suatu sistem dinamik itu sendiri dengan menganalis keadaan disekitar sistem pada keadaan stabil. Diberikan sistem dinamik dua dimensi sebagai berikut: A ′ = @(A, C) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (11) C ′ = (A, C) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (12)

Sebagai contoh untuk menganalisis persamaan berikut S′ = GS + IT ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (16) T ′ = PS + NT ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (17) Bentuk berikut

matriksnya

′ G U S ′V = 8 P T

adalah

sebagai

I S 9 8 9 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (18) N T

matiks linierisasi yang terkait adalah

G I W=8 9 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (19) P N disebut matriks jacobian. 2.3.3 Nilai eigen dan vektor eigen Sebuah vektor yang elemennya tidak ada yang nol disebut vektor eigen V dari sebuah matriks L yang berkaitan dengan nilai eigen λ jika;

8

LV = λV (notasi matriks)………… (20) nilai eigen sangat penting dalam hal analisis sistem dinamik dilihat dari stabilitas titik equilibriumnya. Untuk menentukan nilai eigen harus melalui suatu persamaan karakteristik berikut: G−Y NOX 8 P

karakteristik sebagai kombinasi tiap-tiap kemungkinan nilai eigen seperti pada Gambar 7.3

I 9 = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (21) N−Y

bentuk polinomial matriks diatas adalah

dari

persamaan

(G − Y)(N − Y) − IP = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (22) atau, Y/ − 'Y + ∆= 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (23) ' XMW G+N 8 9=8 9=8 9 ∙∙∙∙∙∙∙∙∙ (24) ∆ det W GN − IP

sebagi suatu fungsi polinomial maka memiliki dua nilai solusi dalam bentuk Y.,/ =

τ±

√' /

√' / − 4∆≥ 0

2

− 4∆

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (25)

Gambar 7. Klasifikasi titik equilibrium.3 •

nilai eigen bernilai real (nyata) jika jika √' / − 4∆< 0 . Pada keadaan ini solusi umum dari sistem linier ini berbentuk U

atau

komplek-konjugat

Titik node, terjadi jika nilai eigen adalah real dan memiliki tanda yang sama. Titik ini stabil ketika nilai eigen keduanya bernilai negatif dan tidak stabil ketika keduanya bernilai positif. Trayektori bersifat konvergen saat stabil dan divergen saat tak stabil Gambar 8.

S(X) V = P. O `.a b. + P/ O `/a b/ ∙∙∙∙∙ (26) T(X)

ketika nilai eigen keduanya bernilai negatif maka akan stabil. Jika sedikitnya satu nilai eigen bernilai positif maka akan tidak stabil. 2.3.4 Klasifiaksi equilibrium Nilai eigen mempengaruhi karakteristik geometri di dekat titik equilibrium. Dari nilai eigen yang diperoleh, didapatkan titik-titik

Gambar 8. Tititk equilibrium node ,λ1=-1, λ2=-3 (stabil), λ1=+1, λ2=+3 (tidak stabil).3 •

Titik saddle, terjadi jika nilai eigen adalah real dan memiliki tanda yang

9

berlawanan. Titik saddle selalu tidak stabil, dikarenakan terdapat nilai eigen yang bernilai positif. Vektor eigen bersifat konvergen menuju nilai eigen negatif dilanjutkan divergen dari nilai eigen positif.

Gambar 9. Saddle equilibrium λ1=+1, λ2=-1.3 •

Titik focus nilai eigen imajiner (kompleks-konjugat). Stabil pada saat nilai eigen memiliki suku real negatif dan tak stabil ketika memiliki suku real yang positif. Suku imajiner dari nilai eigen menentukan frekuensi rotasi dari trayektori yang mengelilingi titik focus.3

Gambar 10.Focus equilibrium λ1=3±i (stabil) atau λ1=+3±i (tidak stabil).3 2.4 Propagasi saraf Jaringan saraf memiliki sifat dapat terangsang (excitability). Dalam kasus ini berarti suatu jaringan saraf memiliki sifat khas pada keadaan istirahat menuju keadaan eksitasi yang menghasilkan suatu potensial aksi pada keadaan ketika saraf mengalami

rangsangan.3 Berbagai macam stimulus dari luar sel saraf sangat mempengaruhi prilaku saraf saat setelah terstimulasi. Saraf yang terkena stimulus, belum tentu akan menghasilkan suatu potensial aksi jika stiumulus tersebut tidak sesuai dengai karakteristik kimia dan fisik dari jaringan saraf tertentu, atau tidak cukup untuk mencapai nilai potensial ambang (threshold potential). Dari sudut pandang sistem dinamik, saraf memiliki sifat excitable atau dapat terangsang karena saraf memiliki bifurkasi keadaan keseimbangan yang dekat dengan keadaan potensial aksi dalam ruang fasenya.3 Tipe dari bifurkasi ini akan menentukan tipe propagasi dalam saraf. Dalam hal ini, akan dijelaskan hubungan antara jenis bfurkasi dengan tipe propagasi saraf yang berkaitan dengan sistem dinamik saraf tersebut. Sebuah sistem dinamik pada saraf dengan dengan sebuah keseimbangan yang stabil, adalah excitable. Ini berarti bahwa jika terdapat suatu rangsangan yang letaknya dekat dengan titik keseimbangan, maka akan meninggalkan titik keseimbangan tersebut (dengan stimulus yang cukup) dan akan membentuk suatu trayektori dan akan kembali lagi ke titik keseimbangan tersebut seperti pada Gambar11.3 Jika stimulus tersebut terjadi secara berkelanjutan, maka akan menimbulkan trayektori yang periodik sehingga membentuk suatu limit cycle, Seperti dapat dilihat pada Gambar 11.

10

•

Gambar 11.Suatu limit cycle. (a)keadaan terangsang excitability.Sebuah trayektori (kotak kecil) terinisiasi dengan titik kestabilan dan kembali ke tiik kestabilan tersebut.(b)sistem tereksitasi dekat dengan bifurkasi akan membenmtuk limit cycle karena titik keseimbangannya belum terinisiasi (selama stimulus diterapkan secara kontinu).3 2.4.1 Bifurkasi Dari bahasan sebelumnya, bifurkasi yang terjadi pada sistem dinamik saraf berkaitan dengan nilai dan vektor eigen. Sebuah sistem dinamik akan stabil apabila semua nilai eigen dari matriks jacobian pada keadaan setimbang memiliki nilai real negatif. Ketika suatu parameter berubah (misal arus eksternal I ) berubah maka nilai eigen akan berubah tanda. Ini dapat menyebabkan suatu bifurkasi saat keseimbangan (bifurcation of the equilibrium).3 Dengan asumsi sistem memiliki dua keadaan stabil (Bistable), Selama keseimbangan dekat dengan limit cycle dan berlangsung secara periodik, Dua kemungkinan bifurkasi yang dapat terjadi adalah sebagai berikut: • Sebuah nilai eigen yang negatif akan berubah menjadi nol. Ini dapat terjadi pada tipe bifurkasi saddle-node dan titik keseimbangan menghilang baik di dalam maupun luar trayektori.

Dua nila eigen kompleks-konjugat dengan bagian real negatif mendekati sumbu imaginer sehingga nilainya akan imaginer. Ini terjadi pada bifurkasi tipe Andronov-Hopf, dan titik keseimbangan akan kehilangan kestabilan tetapi tidak menghilang. Dari kombinasi ketiga titik keseimbangan (node, saddle, focus), dengan perubahan tanda pada nilai eigen, maka akan dihasilkan 4 tipe bifurkasi pada keadaan transisi istirahat (equilibrium) menuju keadaan eksitasi (spiking state). Disajikan pada Gambar 12.

Gambar 12.Transisi terjadi dari keadaan istirahat menuju keadaan eksitasi yang terjadi melalui bifurkasi keseimbangan (tanda panah).Bifurkasi saddle-node pada in vitropyramidal neuron of rat’s primary visual cortex.BifurkasiAndronov-Hopf pada in vitro brainstem mesencephalic V neuron.3 Menurut sebuah teori yang berkaitan dengan sistem dinamik pada saraf yang menghubungkan sebuah sistem tereksitasi adalah merupakan bagian dari suatu bifurkasi dari keadaan

11

dinamik (spiking state) yang berulang (oscillatory dynamics) menuju kedaan istirahat (resting) akan menghasilkan suatu bifurkasi pada keadaan eksitasi (bifurcation of a limit cycle). Bifurkasi ini terjadi ketika sebuah trayektori terinisiasi diluar titik keseimbangan dan memiliki kemampuan untuk eksitasi secara periodik. Ini akan menghasilkan sistem dinamik yang menampilkan suatu limit cycle (Gambar 13).

Gambar 14. Sistem dinamik tereksitasi dapat memiliki bifurkasi yang berbeda baik melalui keadaan bistable atau tidak.3

Gambar 13.Limit cycle empat bifurkasi pada keadaan istirahat.3 Dalam hal ini, suatu sistem dinamik dapat beralih dari keadaan istirahat (resting state) menuju keadaan eksitasi (spiking state) baik melalui keadaan bistable maupun tidak (monostable). Secara singkat dapat dilihat pada Gambar 14.

2.4.2 Klasifikasi propagasi saraf Mekanisme bifurkasi pada saraf tereksitasi pertama kali dipelajari oleh Hodgkin (1948). Ia menerapkan berbagai variasi nilai arus pada membran tereksitasi sel saraf otak tikus dan mencatat hasil yang didapatkan. Dalam penelitiannya menunjukan ketika arus yang diterapkan lemah, saraf target tidak menunjukan reaksi apapun. Ketika arus yang diterapkan kuat, maka terjadi eksitasi pada membran sehingga menyebabkan potensial aksi. Berdasarkan besar atau kecilnya nilai rata-rata arus yang diterapkan pada membran untuk terjadinya suatu potensial aksi, secara garis besar

12

Hodgkin (1948) mengklasifikasikan eksitasi saraf menjadi dua kelas.24 • Eksitasi saraf tipe 1 (class 1 neural excitability). Suatu potensial aksi dapat dihasilkan dengan frekuensi acak rendah dan bergantung pada arus yang diterapkan. • Eksitasi saraf tipe 2 (class 2 neural excitability). Suatu potensial aksi dapat dihasilkan dalam suatu pita frekuensi tertentu dan relatif tidak bergantung terhadap perubahan kekuatan arus yang diterapkan.

Gambar 15. Propagasi eksitasi saraf (a) tipe 1 dan (b) tipe 2 dengan arus terapan DC konstan pada sela saraf tikus.3 Tipe eksitasi 1 dan 2 ini lah yang akan dibahas pada penelitian ini dengan perlakuan variasi nilai arus terapan sebagai suatu konstanta, arus DC atau pun AC bergantung waktu. Selain itu, dibahas pula untuk analisis dinamik dan model saraf terkopel untuk kasus yang lebih kompleks. 2.5 Sinkronisasi Otak manusia merupakan sebuah jaringan kompleks yang memiliki sekitar 1011 sel saraf yang saling terhubung oleh

1014 sampai 1015 penghubung 10 (conector). ini menunjukan bahwa fenomena sinkronisasi antar sel saraf dalam proses pengolahan informasi dalam bentuk impuls adalah suatu hal yang pasti terjadi dalam sistem dinamik pada saraf.11 Sistem kompleks pada saraf merupakan suatu set kumpulan unit saraf yang terhubung dan saling terosilasi yang dapat dimodelkan sebagai sutau unit sel yang berpasangan (coupled oscillator) . coupled oscillator ini saling tersinkronisasi dalam proses pengolahan informasi pada saraf.12 Salah satu cara untuk menganalisis dan memahami fenomena ini adalah dengan mengasumsikan sistem kompleks ini menjadi suatu unit kecil osilator, kemudian menganalisis mode ini pada dimensi yang lebih kecil dari sudut pandang fisisnya saja. Hoppensteadt dan Izhikevich (1997) telah melakukan beberapa analisis kualitatif pada sistem kompleks ini dengan memodelkan suatu solusi periodik dan sinkronisasinya pada dua sel saraf terkopel.3 Analisis yang dilakukan oleh Izhikevich dan Hoppensteadt adalah dengan mereduksi model satu saraf menjadi suatu model fase (phase model) yang diturunkan dari sistem dinamik limit cycle pada spiking neuron.. Hasil dari analisis model fase dapat digunakan untuk mengetahui karakteristik dari osilator tekopel pada model saraf kompleks.12 Model fase digambarkan oleh suatu variabel fase (ϑ). Θ didefinisikan sebagai fase setelah spiking neuron. Dua osilator terkopel dapat menempati tiga keadaan sinkronisasi yaitu in-phase, antiphase, atau out-of-phase ketika perbedaan fasenya ϑ2- ϑ1 masing-masing sebanding dengan 0, setengan periode, atau variasi nilai lainnya. (Lihat Gambar 16).3

13

′

= @(b ) + q g f

g

h.

jbg m ∙∙∙∙∙∙∙∙∙∙ (28)

K = 1,2,3, … . . , os = 1,2,3, … … o

(a) (b) (c) Gambar. 16. Tiga kondisi sinkronisasi (a) sefase (b) berlawanan fase (c) berbeda fase.3 Secara umum, osilator saraf terkopel dapat di modelkan seperti persamaan berikut: A ′ = @ (A ) + e f gh.

g

(A )ijX − Xg∗

− l g m ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (27)

Dengan xi € [0,1] merupakan potensial membran ke-I pada saraf. Fungsi fi merupakan persamaan pada sistem dinamik.. ketika xi mencapai 1, maka saraf ke-I akan spiking dan xi akan kembali ke 0. Xg∗ merupakan waktu yang dibutuhkan saat spiking hingga reset. Parameter dimensional ε menggambarkan kekuatan dari koneksi antar unit saraf. Fungsi gij mewakili parameter yang berkaitan dengan efek dari pasangan unit saraf ke-j saat spiking terhadap unit saraf ke-I saat spiking. ε gij(xi) merupakan increament bagi xi setelah lij ≥ 0. Increamnet ini dihasilkan dari suatu fungsi delta dirac δ dengan δ(t) =0 untuk semua X ≠ 0, i(0) = ∞NGo p i = 1 dengan mengasumsikan fi dan gij adalah kontinu.13 Bentuk sederhana dari persamaan (27) adalah model n saraf dengan pengaruh keterhubungan antar unit saraf dan kekuatannya.

V merupakan potensial membran, f(vi) merupakan fungsi potensial membran kei seperti pada persamaan (21). ε merupakan parameter yang bertanggung jawab atas kekuatan kopling. Sedangkan gij merupakan sutau fungsi kopling antar saraf. Fungsi kopling ini memiliki berbagai macam bentuk bergantung terhadap acuan apa model fungsi tersebut dibangun. Pada model ini, diasumsikan bahwa saraf terhubung satu dengan lainnya secara sinaptik dengan fungsi kopling sebagai fungsi sigmoid berikut.4 g jbg m

=

1

1 + exp 8σjbg − tg m9

∙∙∙∙ (29)

Dengan parameter σ berperan dalam pengaturan laju kopling dan θj merupakan potensial ambang tiap sel saraf ke-j. Contoh hasil simulasi yang telah dilakukan oleh Izhikevich dengan menggunakan model integrate and fire untuk dua saraf terkopel disajikan pada Gambar 17.13

Gambar 17. Dua saraf terkopel dengan peta frekuensi impuls (a)berbeda fase (b) sefase (c) berlawanan fase.31

14

Pada model integrate and fire terkopel memiliki dua kemungkinan keadaan sinkronisasi yaitu sefase (atas) dan berbeda fase (bawah). Anaisis sinkronisasi ini dapat diterapkan pada berbagai model seperti Hodgkin-Huxley, Fitzhugh-Nagumo, dan Morris-Lecar karena merupakan suatu set persamaan sistem dinamik dimensional.14

III. METODE PENELITIAN 3.1 Waktu dan Tempat Penelitian Penelitianinidilakukan di Laboratorium Fisika Teori dan Komputasi Departemen Fisika Institut Pertanian Bogor (IPB) dan kediaman peneliti. Waktu penelitian dilakukan selama 5 bulan dari bulan 8 Desember 2010 hingga 8 Mei 2011. 3.2 Peralatan Peralatan yang digunakan dalam penelitian ini adalah sebuah PC dengan prosesor Intel Core 2 Quad (2.3 GHz), 4 GB RAM (Lab. Teoridan Komputasi) dan PC dengan prosesor Intel Core 2 Duo (2.9 GHz), 2 GB RAM (Kediaman peneliti). Perangkat lunak yang digunakan dalam penelitian ini adalah MS. Office 2007 dan MATLAB R2009b. Selain PC, peralatan yang digunakan berupa bahan rujukan dan pustaka yang peneliti dapat kandari internet, perpustakaan dan buku catatan. 3.3 Metode Penelitian 3.3.1 Studi pustaka Penelitian ini diawali dengan studi pustaka yang meliputi pencarian pustaka mengenai mekanisme transport rmembran pada saraf, penjalaran impuls saraf, hingga pustaka mengenai model saraf dalam persamaan matematis. Pemecahan model saraf dalam bentuk persamaan differensial(PD) dilakukan dengan analisis numerik menggunakan metode Rungge-Kutta4 (RK-4) dengan simulasi MATLAB pada berbagai variasi parameter. Sedangkan untuk karakteristik dinamiknya dipecahkan dengan cara analisis sistem dinamik dengan

menerapkan nilai eigen dari matriks jacobian. 3.3.2 Analisis numerik solusi propagasi saraf Morris-Lecar dengan metode RK-4 Suatu fungsi dalam differensiasi numerik diperlukan dalam perkembangan penggunaan algoritma pemrograman untuk memecahkan suatu nilai batas untuk PD biasa (PDB) dan PD parsial (PDP). Contoh yang biasa digunakan pada differensiasi numerik yaitu pada penggunaan pendekatan numerik yang dapat dibandingkan dengan solusi eksak. Salah satu metode pemecahan PD secara numerik yaitu dengan menggunakan metode RK-4. RK-4 merupakan metode pendekatan analisis PD numerik yang memiliki tingkat ketelitian yang tinggi. Dengan nilai increament yang sangat kecil dan berulang-ulang kemudian diambil nilai rata-ratanya, dapat memberikan tingkat pendekatan dengan hasil eksak yang akurat. Model saraf Morris-Lecar (1981) terdiri dari dua PDB dengan duavariabel dimensi (V) dan (W) . selanjutnya model PDB ini dianalisis dengan menggunakan RK-4 untuk menampilkan solusi secara numerik. Solusi numerik yang akan dianalisis dilakukan pada arus terapan tetap, DC bergantung waktu, dan AC bergantung waktu. 3.3.3 Analisis sistem dinamik model saraf Morris-lecar Analisis sistem dinamik ini dilakukan dengan mencari nilai eigen, dan nilai akar-akar nol pada model PDB ML. Selanjutnya dibangun suatu matriks Jacobi untuk mengetahui jenis titik kritis pada sistem yang dapat menetukan jenis bifurkasi pada model PDB dengan analisis pada ruang fase pada sistem. Analisis sistem dinamik ini dilakukan pada arus terapan tetap, DC bergantung waktu, dan AC bergantung waktu. 3.3.4 Sinkronisasi model saraf MorrisLecar terkopel Proses ini merupakan simulasi terakhir dari penelitian ini dengan asumsi bahwa sistem saraf kompleks dapat dimodelkan oleh suatu sistem saraf

15

terkopel. Analisis yang dilakukan pada sistem terkopel ini dilakukan hanya pada model dengan arus AC bergantung waktu saja. Pada sistem terkopel ini akan dibahas propagasi sistem kompleks saat terisolasi (tidak terkopel), terkopel, dan sinkronisasi dengan variasi fase propagasi yang berbeda dengan melibatkan kekuatan kopel antar saraf. Hasil yang didapat dalam analisis ini ditampilkan dengan menggunakan MATLAB berupa propagasi sistem banyak saraf (n=2,3,4) terkopel.

IV. HASIL DAN PEMBAHASAN Selama beberapa dekade terakhir ini, penelitian mengenai jaringan saraf tiruan (JST) berkembang seiring dengan kemajuan berbagai teknologi perangkat lunak dalam hal analisis JST tersebut. Dalam proses interpretasi JST, berbagai model telah dipublikasikan oleh para peneliti untuk memvisualisasikan bagaimana mekanisme propagasi pada jaringan saraf dalam bentuk action potential (AP). Salah satu model yang telah berhasil memvisualisasikan mekanisme AP pada jaringan saraf adalah model Morris-Lecar (1948) yang merupakan sistem pesamaan differensial biasa (PDB) terhadap waktu dengan dua variabel dimensional utama yaitu V dan W. Dengan meninjau kembali persamaan (2) dan (3), model saraf Morris-Lecar (ML) merupakan model yang diaplikasikan untuk suatu sistem jaringan saraf yang memiliki sensitifitas terhadap tegangan listrik akibat adanya konduktansi pada membran sel saraf.15 Model ini memiliki dua variabel dimensional utama yaitu V dan W yang masing-masing mewakili potensial membran saraf dan suatu recovery variable yang berhubungan dengan normalisasi konduktansi ion K+ dalam peristiwa depolarisasi. Fungsi ini dibangun berdasarkan asumsi bahwa nilainya sebanding dengan nilai instan dari kemungkinan saluran ion tersebut berada pada keadaan terbuka. Iapp merupakan variabel yang bertanggung

jawab atas adanya rangsangan dari luar berupa arus listrik yang diterapkan pada sel saraf. C merupakan parameter kapasitansi total dari membran saraf. Parameter VCa, VK,, dan Vl mewakili potensial kesetimbngan dari ion Ca2+, K+, dan faktor koreksi dari arus kebocoran (Leakage Current). Sedangkan gCa, gK, merupakan konduktansi dan gl, maksimum yang bertanggung jawab atas arus ionik yang terjadi pada sel saraf. Fungsi M∞(V) bergantung pada nilai potensial membran merupakan sutau fungsi yang berkaitan dengan peluang terbukanya saluran Ca2+ dapat dilihat pada persamaan (4). Persamaan (5) menggambarkan proses pemulihan yang dilakukan oleh saluran protein yang bertransformasi dengan membran saraf diantara keadaan terkonduksi ion-ion atau tidak. Pada persamaan kedua ini terdapat dua buah fungsi kemungkinan W∞dan τ∞ yang masing masing merupakan fungsi kemungkinan terbukanya saluran K+ dan suatu fungsi skala waktu yang berkaitan dengan proses pemulihan (depolarisasi). Pada persmaan (8), parameter ø merupakan skala waktu proses pemulihan. Nilai ø dapat divariasikan untuk berbagai sel yang berbeda-beda dan sangat sensitif terhadap suhu lingkungan membran. Parameter V1, dan V3 merupakan suatu nilai tengah saat arus ionik Ca2+ dan K+ada pada keadaan setengah teraktivasi (half activated), V2 merupakan sebuah konstanta potensial yang bertanggung jawab kepada loncatan potensial saat aktivasi, sedangkan V4 adalah faktor kemiringan laju aktivasi ion K+.16 Secara keseluruhan, saat saraf menerima rangsangan dari luar ,maka akan terjadi suatu potensial aksi karena mekanisme elektrik yang menyebabkan perubahan beda potensial, arus, konduktansi, dan kapasitansi pada membran dalam proses penjalaran impuls tersebut. 4.1 Solusi Numerik Propagasi Saraf dengan Metode RK-4 Untuk menyelesaikan PDB diatas digunakan pendekatan secara numerik

16

Dalam pendekatan secara numerik, solusi yang akan dibangun merupakan hasil iterasi PDB dengan anggapan bahwa nilai V dan W akan berubah terhadap selang waktu dt. Sehingga dalam hal ini variabel dt merupakan suatu parameter iterasi pada suatu pendekatan numerik atau sering disebut sebagai increament. Persamaan (30) dan (31) dapat disederhanakan penulisannya menjadi suatu fungsi f(v,w) dan g(v,w).dengan membuat ruas kiri kedua persamaan masing-masing hanya terdiri dari parameter dV dan dW, maka persamaan sebelumnya akan menjadi persamaan (32) dan (33), N = @(b, T)NX ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (32) N! = (b, T)NX ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (33)

dengan nilai f(v,w) dan g(v,w) masingmasing: @(b, T) = (− _ G _∞ ( )( − _ G ) − _w !( − _w ) − _W ( − _W ) + $_Gxx)/ ∙∙∙∙∙∙ ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (34) !∞ ( ) − ! (b, T) = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (35) '( ( )

Bentuk persamaan (32) dan (33) ini dianalisis secara numerik (Lampiran 3) dengan menggunakan metode RK-4. 4.1.1 Solusi numerik dengan arus terapan DC tetap

Dengan menggunakan perangkat lunak MATLAB, didapatkan hasil analisis numerik pada model ML yang disajikan pada Gambar 18. Dengan nilai parameter yang terkait adalah C=20 µF/cm2, gK=8 ms/cm2, gl=2 ms/cm2gCa=4 ms/cm2 , ø=1/15 s-1, VCa= 120 mV, VK=80 mV, Vl= -60 mV, V1=-1.2 mV, V2=18 mV, V4=17.4 mV ,V3=12 mV. dan Iapp= 50 µA. 40

30

20

m e m b ra n e v o lta g e v (m V )

dengan menggunakan metode RunggeKutta orde-4 (RK-4). V’ merupakan nilai perubahan potensial membran terhadap waktu yaitu dV/dt sedangkan W’ merupakan laju proses depolarisasi pada membran dW/dt sehingga persamaan (2) dan (3) menjadi. N ) =− ∞ ( )( − NX − !( − ) − " ( − " ) + $ %% ∙∙∙∙∙ ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (30) N! !∞ ( ) − ! = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (31) NX '( ( )

10

0

-10

-20

-30

-40

-50

0

200

400

600 time t (ms)

800

1000

1200

Gambar 18. Aktivitas listrik (action potential) model saraf Morris-Lecar tipe 1 Program dengan metode RK-4 dapat dilihat selengkapnya pada Lampiran 2. Pada bab 2, telah dijelaskan bahwa mekanisme propagasi saraf memiliki berbagai macam bentuk dinamik (neural properties). Dalam hal ini, antara sel satu dengan yang lain memiliki karakteristik spesifik saat menerima rangsangan dari luar. Baik ditinjau dari kecepatan responnya, besar kecil rangsangan (applied current) , nilai resting potential (RP), maupun sifat dinamik dalam propagasinya. Semua kombinasi ini menghasilkan suatu mekanisme dinamik yang bervariasi dalam suatu propagasi saraf. Bentuk propagasi yang dibahas dalam penelitian ini seperti yang telah di klasifikasikan oleh Hodgkin (1948) dilihat dari segi rata-rata frekuensi arus yang diterapkan pada sel untuk suatu peristiwa eksitasi adalah Eksitasi Saraf Tipe 1 (class 1) dan Eksitasi Saraf Tipe 2 (class 2). Gambar 1.merupakan bentuk propagasi class 1 dengan nilai arus Iapp merupakan arus DC dengan nilai yang konstan. Dengan menggantikan nilai

17

parameter V3 menjadi 2 mV dan Iapp= 55 µA maka didapatkan bentuk propagasi class 2 seperti pada Gambar19.

30 20 Class 2

(b)

10 0

40 -10

Class 2 Excitability

30

-20

membran voltage v (mV)

20

-30

10

-40

0

-50

0

100

200

300

400

500

600

-10 -20 -30 -40 -50

0

200

400

600 time t (ms)

800

1000

1200

Gambar 19. Bentuk propagasi saraf tipe 2. Hasil simulasi tidak menunjukan adanya perbedaan antara Tipe 1 dan 2. Kedua tipe propagasi tersebut sebenarnya memiliki perbedaan dalam hal sistem dinamiknya. Perbedaan nilai titik keseimbangan dan jenis bifurkasi sangat jelas terlihat pada suatu bidang fase pada tipe 1 dan 2. Pembahasan lebih lengkapnya, akan dijelaskan pada sub bab berikutnya. Berdasarkan hasil simulasi, pada kedua tipe propagasi memiliki nilai minimum Iapp untuk melakukan eksitasi secara periodik (Gambar 18 dan 19). Nilai minmum untuk tipe 1 dan 2 masing-masing adalah 40 mA dan 50 mA. nilai ini merupakan nilai minimum agar suatu potensial aksi dapat menjalar secara periodik. Jika nilai Iapp≤ Imin, maka sel saraf tersebut tidak cukup kuat untuk mengirimkan sinyal, atau dalam arti lain hanya mampu melakukan sekali eksitasi kemudian akan kembali ke keadaan istirahat. 40 30

Class 1

20 10

(a)

0 -10 -20 -30 -40 -50

0

100

200

300

400

500

600

Gambar 20. Nilai Iapp pada (a) tipe 1 dan (b) tipe 2 masing-masing 40 µA dan 50 µA. Kedua bentuk propagasi tidak dapat terjadi secara periodik. 4.1.2 Solusi numerik dengan arus terapan DC bergantung waktu Nilai arus Iapp atau arus yang diterapkan pada sel saraf sangat mempengaruhi bentuk propagasinya. Pada sub bab sebelumnya, telah dibahas bentuk propagasi saraf pada tipe 1 dan 2 dengan nilai arus terapan adalah konstan, yaitu masing-masing 50 µA dan 55 µA untuk tipe 1 dan 2. Dengan nilai tersebut, saraf dapat menjalar secara periodik. Jika arus Iapp pada sel saraf tidak bernilai tetap, atau nilainya berubah terhadap waktu, maka bentuk propagasi dan sistem dinamiknya berubah. Dalam penelitian ini dimodelkan suatu persamaan yang merupakan fungsi arus terapan Iapp terhadap waktu I(t) sebagai berikut: $(X) = $z { |X + $

a

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (36)

fungsi arus I(t) pada persamaan (36) dimodelkan sebagai suatau fungsi linier yang berbanding lurus dengan waktu. Ini berarti bahwa nilai arus terapan pada sel saraf akan berubah dengan bertambahnya waktu. Parameter Imax merupakan nilai penambahan (gradien) arus maksimum tiap detik, sedangkan α merupakan nilai koefisien penambahan yang bertanggung jawab atas besar kecil laju perubahan arusnya. Dengan mensubstitusikan persamaan (36) ke persamaan (34) dengan menggantikan parameter Iapp dengan I(t), persamaan (34) menjadi persamaan (37) sebagai berikut:

18 @(b, T) = (− _ G _∞ ( )( − _ G ) − _w !( − _w ) − _W ( − _W ) + $(X))/ (37)

Propagasi saraf tipe 1 dan 2 ini memiliki karaktersitik masing-masing dalam merespon rangsangan dari luar. Dengan mengubah nilai Iapp menjadi suatu nilai yang bergantung dengan waktu, Nilai parameter kedua tipe berbeda. Selain I(t), nilai V3 padakedua tipe berbeda yaitu 12 mV dan 2 mV untuk tipe 1 dan 2. Perbedaan nilai ini pada kedua tipe saraf tersebut menampilkan bentuk propagasi yang berbeda. Berdasarkan Gambar 21., tipe 1 mulai melakukakn eksitasi pada saat t≈800 ms (spike state) yaitu pada saat nilai I≈130 µA. Saat nilai I sangat besar (I≈350 µA) potensial aksi mulai menghilang (t≈2050 ms). Sedangkan untuk tipe 2 (Gambar 22.) saraf mulai tereksitasi saat t≈350 ms dengan nilai I≈60 µA dan saat t≈1600 ms (I≈260 µA) propagasi berada pada keadaan istirahat. Kondisi ini berkaitan dengan karaktersitk saraf. Sebagai suatu komponen biologi fungsional, sel saraf memiliki karakteristik spesifik dalam merespon rangsangan dari luar. Secara fisis, sel-sel saraf pada tubuh cenderung sensitif terhadap adanya rangsangan dari luar berupa adanya arus yang diterapkan. ketika nilai arus yang diterapkan tidak cukup untuk melakukan depolarisasi maka tidak akan terjadi suatu potensial aksi. Ketika mulai mencapai potensial ambang, maka akan terjadi suatu potensial aksi. Jika nilai arus yang diterapkan melebihi ambang batas saraf, atau diluar interval saraf untuk menghasilkan suatu potensial aksi, maka tidak akan terjadi propagasi pada saraf.18

persamaan (37) kemudian disubstitusikan kembali ke persamaan (32), kemudian dengan menggunakan MATLAB didapatkan solusi numerik seperti pada Gambar 21. (class 1) dan Gambar 22. (class2). Pulse of Class 1 Current Time Dependent

40

m e m b ra n e v o lta g e (m V )

20

0 rest state

spike state

-20

-40

-60

0

500

1000

1500

2000

2500

1500

2000

2500

time (ms)

a p p lie d c u rr e n t (m ic ro A m p e re )

400

200

0

0

500

1000 time (ms)

Gambar 21. Propagasi saraf tipe 1 dengan arus I(t). parameter untuk propagasi tipe 1 adalah Imax= 5 µA, Iinit= 0, dan α=0.011 s-1, sedangkan untuk tipe 2 adalah Imax= 10 µA, Iinit=0, dan α=0.016 s-1. Pulse of Class 2 current time dependent 40

m e m b ra n e v o lta g e (m V )

20

0 spike state

rest state

-20

30

-40

20

10

-60

spike state

0

-80

0

200

400

600

800

1000 time (ms)

1200

1400

1600

1800

-10

-20

400 a p p l ie d c u r r e n t (m ik ro A m p e re )

rest state

2000

Class 2

-30

200 -40

0

0

200

400

600

800

1000 time (ms)

1200

1400

1600

1800

2000

-50

(a)

36 spikes/1200 ms -60

Gambar 22. Propagasi saraf tipe 2 dengan arus I(t).

-70

0

500

1000 time (ms)

1500

2000

19

(b)

0

res t s tate

-10

s pike s tate

-20

-30

-40

-50

28 s pikes/1200 ms -60

0

500

1000 1500 time (ms )

2000

2500

Gambar 23. Frekuensi Frekuensi propagasi (spike/second) pada (a) tipe 1 dan (b) tipe 2 Pada model ini, kedua tipe saraf tersebut memiliki nilai resting potential yang hampir sama yaitu sekitar -60 mV. Bentuk propagasi saraf tipe 1 dan 2 merupakan tipe eksitasi saraf utama yang digolongkan berdasarkan besar atau kecilnya nilai rata-rata arus yang diterapkan pada membran untuk terjadinya suatu potensial aksi. Hodgkin (1948) menklasifikasikan bahwa propagasi tipe 1 dapat dihasilkan dengan frekuensi eksitasi yang rendah dan bergantung pada besar arus yang diterapkan. Sedangkan untuk tipe 2 dapat terjadi hanya pada pita frekuensi eksitasi tertentu dan tidak bergantung oleh besar arus yang diterapkan. Berdasarkan hasil yang ditampilkan pada Gambar 23., dapat dilihat bahwa frekuensi eksitasi pada tipe 2 (36 spikes/1200 ms) lebih besar dari tipe 1 (28 spikes/1200 ms). Berdasarkan hasil eksperimen Hodgkin (1848) dan penelitian lebih lanjut oleh E. M. izhikevich (2003), menunjukan bahwa perbedaan kualitatif antara tipe 1 dan 2 ditandai oleh nilai arus yang diterapkan pada sel. Arus terapan akan kontinu dan menuju stabil dalam menghasilkan suatu potensial aksi untuk tipe 1, sedangkan tipe 2 memiliki nilai rentang arus tertentu untuk menghasilkan suatu potensial aksi. Jika di luar pita ini, maka tidak dapat dihasilkan suatu potensial aksi. Agar lebih memahami teori pita frekuensi pada eksitasi tipe 1 dan 2, akan ditinjau kembali nilai I(t). Nilai Iapp pada

40 30

Class 1

20

(a)

10

m em brane v oltage (m V)

10

0 -10 Resting State

-20 -30 -40 -50 Periodic Spike -60

0

200

400

600

800

1000

1200

1400

1600

1800

1000

1200

1400

1600

1800

time (ms) 100 applied curent (m icroA m pere)

Class 1 20

50 0 -50 -100

0

200

400

600

800 time (ms)

60 40

(b)

no spike m em bran e v oltag e (m V )

30

20

0 no spike

-20

Class 2 -40

-60 Periodic Spike -80

0

200

400

600

800

1000

1200

1400

1600

1800

1000

1200

1400

1600

1800

time (ms) 100 ap plied c urrent (m ic ro A m p ere)

40

model sebelumnya memiliki gradien yang positif bahwa nilai arus akan semakin meningkat dengan bertambahnya waktu. Parameter yang bertanggung jawab dalam hal ini adalah α yang bertanda positif (+). Dengan mengubah tanda pada parameter α menjadi negatif (-), maka gradien fungsi akan negatif sehingga menyebabkan fungsi arus terapan akan terus berkurang dengan bertambahnya waktu. Dengan menggunakan nilai parameter sebelumnya dan mengubah nilai Iinit pada tipe 1 dan 2 masing-masing bernilai 100 µA dan 280 µA, maka didapatkan bentuk propagasi seperti pada Gambar 24.

0 -100 -200

0

200

400

600

800 time (ms)

Gambar 24. Propagasi (a) tipe 1 dan (b) tipe 2 dengan gradient I(t) negatif Teori mengenai propagasi tipe 1 dan 2 dapat dijelaskan dengan melihat hasil yang didapatkan pada Gambar24. Pada tipe 1, proses eksitasi periodik terus terjadi bersamaan dengan perubahan nilai arus Iapp, hingga pada nilai Iapp tertentu saraf tidak cukup energi untuk melakukan eksitasi karena nilai Iapp yang terus berkurang. Sedangkan pada tipe 2, pita frekuensi eksitasi terlihat dengan jelas. Eksitasi saraf periodik hanya terjadi pada pita frekuensi tertentu yaitu pada selang sekitar 500-1500 ms,

20

$(X) = $z { sin (~X) + $

a

∙∙∙∙∙∙∙∙∙∙∙ (38)

Dengan mengganti fungsi I(t) pada persamaan (37) dengan persamaan (38), maka arus terapan pada model akan berupa arus AC yang nilainya menunjukan suatu hubungan sinusoidal terhadap waktu. Parameter Imax dan Iinit memiliki arti fisis yang sama dengan fungsi arus DC bergantung waktu pada sub bab sebelumnya, sedangkan parameter yang berbeda adalah ω yang merupakan nilai frekuensi masukan pada sinyal arus AC yang diterapkan pada model. Dengan memasukan nilai Imax , Iinit dan ω pada tipe 1 dan 2, maka dihasilkan suatu propagasi saraf seperti Gambar 25.

m e m b r a n e V o lt a g e ( m V )

30 20 10

0

-10 -20 -30 -40 -50

0

500

1000

1500

2000

2500

3000

2000

2500

3000

a p p lie d c u r r e n t ( A C )

time (ms)

10

0

-10

0

500

1000

1500

time (ms)

(a)

Class 2 excitability with applied AC current

40 30 20

m e m b r a n e v o lt a g e ( m V )

4.1.3 Solusi numerik dengan arus terapan AC bergantung waktu Nilai parameter Iapp dapat divariasikan bedasarkan karakteristik dari tiap-tiap sel pada jaringan saraf. Pada sub bab ini, akan digunakan suatu nilai arus terapan yang bergantung terhadap waktu I(t) dan nilainya selalu berubah. Parameter yang digunakan ini adalah nilai Iapp dengan fungsi masukan berupa nilai arus AC (alternating current) yang dapat dilihat pada persamaan (38).

Class 1 excitability with applied AC current 40

10

0

-10 -20 -30 -40 -50 -60

0

500

1000

1500

2000

2500

3000

2000

2500

3000

time (ms) a p p lie d c u r r e n t ( ( m ik c r o A m p e r e )

dengan nilai Iapp sekitar 50 µA hingga 150 µA. Kedua keadaan diatas, yaitu ketika kedua tipe diberi arus terapan yang berubah terhadap waktu (baik bertambah maupun berkurang) yang artinya bahwa kedua tipe propagasi tersebut memiliki perbedaan dalam sistem dinamiknya. Hal yang harus digaris bawahi adalah, parameter yang diubah pada pendekatan numerik ini hanya parameter-parameter yang berkaitan dengan nilai arus terapan. Jika parameter-parameter diluar arus terapan divariasikan nilainya, maka akan menghasilkan pola propagasi dan sistem dinamik yang berbeda.

20

0

-20

0

500

1000

1500

time (ms)

(b) Gambar 25. Propagasi saraf dengan fungsi arus terapan AC.(a) tipe 1.(b) tipe 2. nilai paramer untuk tipe 1 adalah Imax = 8 mV, Iinit =50 mA dan ω = 0.011 s-1, Sedangkan untuk tipe 2 adalah Imax =10 mV, Iinit = 55 mA dan ω = 0.0016 s-1. Pengaruh adanya masukan arus AC pada kedua tipe propagasi menyebakan perubahan mekanisme sistem dinamik pada masing-masing tipe propagasi. Tipe 1 merupakan propagasi saraf yang dapat mengalami eksitasi saat arus yang diterapkan berada pada frekuensi yang rendah sedangkan pada tipe 2 relatif sedikit lebih tinggi untuk mengalami eksitasi dan memiliki pita frekuensi eksitasi tertentu. Jika dilihat

21

hasil pada Gambar 25., saat nilai arus definit positif, pada tipe 1 maupun 2 mengalami eksitasi. Perbedaan pada kedua tipe propagasi ini terletak pada saat nilai arus masukan bernilai negatif. Pada tipe 1, meskipun nilai arus masukan memasuki negatif, eksitasi masih dapat terjadi tetapi mengalami penurunan frekuensi eksitasi (spike frequence) dibandingkan saat nilai arus adalah positif. Hal yang berbeda terjadi pada tipe 2. Saat nilai arus negatif, pada tipe 2 tidak terjadi eksitasi sama sekali. Ini berkaitan dengan karakteristik dari propagasi tipe 2, karena pada tipe ini saraf cenderung harus diterapkan oleh nilai arus yang lebih tinggi dengan pita frekuensi eksitasi yang lebih sempit (spesifik).17 Agar lebih memahami fenomena ini, pada tiap tipe 1 dan 2 diperlakukan suatu variasi nilai ω. Nilai ω menunjukkan besar kecilnya frekuensi arus listrik masukan AC pada saraf. Nilai variasi ω dapat dilihat pada Gambar 26. 50

omega

ω

0

0.011 -50

0

200

400

600

800

1000

1200

1400

1600

1800

0

200

400

600

800

1000

1200

1400

1600

1800

0

200

400

600

800

1000

1200

1400

1600

1800

200

400

600

800

1000

1200

1400

1600

1800

50

0.051

0

-50 50

0.101

0

-50

50

0.201

0

-50

0

Class 1 Excitability 50

omega 0

0.016

ω

-50

0

200

400

600

800

1000

1200

1400

1600

1800

0

200

400

600

800

1000

1200

1400

1600

1800

0 50

200

400

600

800

1000

1200

1400

1600

1800

50

0.056

0 -50

-100 50

0.106

0 -50

-100

0.206

0 -50 -100

0

200

400

600

800

1000

1200

1400

1600

Class 2 Excitability

Gambar 26. Variasi nilai ω terhadap bentuk propagasi saraf

1800

Berdasarkan hasil simulasi pada Gambar 26., pada propagasi tipe 1, semakin besar nilai ω, perubahan frekuensi spike tidak terlalu besar namun terdapat perubahan fase propagasi menuju stabil. Sedangkan pada tipe 2, perubahan nilai ω yang semakin besar, sangat terlihat perubahan yang signifikan. Pada nilai ω=0.016, tipe 2 melakukan burst, saat nilainya dinaikan menjadi 0.056, propagasi burst menghilang dan menjadi suatu tonic spiking. Saat nilai ω dinaikan lagi menjadi 0.106, peristiwa burst kembali muncul dan saat ω bernilai 0.206 propagasi kembali stabil (regular spiking). Dapat disimpulkan bahwa pada tipe 1, kenaikan nilai ω cenderung tidak mengubah bentuk propagasi saraf (neural properties) hanya mengubah keteraturan propagasi saraf dilihat dari fase propagasi tiap eksitasi (spike) hingga mencapai kestabilan. Sedangkan pada tipe 2, perubahan (kenaikan) nilai ω dapat mengubah bentuk propagasi saraf baik itu berupa spike atau burst secara berulang. 4.2 Analisis Sistem Dinamik Propagasi Saraf Langkah terakhir dari analisis kualitatif suatu sistem dinamik adalah analisis bifurkasi. Suatu sistem dinamik dikatakan mengalami bifurkasi alamiah ketika ruang fasenya memiliki karakteristik perubahan secara kualitatif.3 Perubahan secara kualitatif adalah perubahan karakteristik sistem dinamik saat ada atau tidak ada dalam keadaan dinamik. Suatu sel saraf berada pada keadaan ada atau tidak dinamik bergantung pada kondisi awal dan parameter alamiah yang berkaitan dengan saraf tersebut. Dalam hal ini yang paling terlihat jelas adalah parameter potensial membran. Bifurkasi merupakan proses perubahan titik keseimbangan (equilibrium) baik jenis maupun jumlah akibat adanya perubahan parameter yang terkandung pada suatu persamaan.7 Dalam hal ini parameter dan persamaan yang dimaksud terangkum dalam sutau model saraf. Model yang digunakan adalah model ML dengan parameter

22

)( − ) − !( − ) − " ( − " ) + $ %% = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (41) !∞ ( ) − ! = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (42) '( ( )

−

∞(

ruas kiri pada masing-masing persamaan dimodifikasi sehingga hanya mengandung parameter w saja sehingga persamaan (41) dan (42) menjadi. !( ) = (− _ G _∞ ( )( − _ G ) − _W ( − _W ) + $_Gxx)/( _w ( − _w ) ) ∙∙∙∙∙∙∙∙∙∙∙∙∙ (43) !( ) = !∞ ( ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (44)

Persamaan (43) merupakan grafik garis nol (nullcline) saat nilai dV/dt=0 sedangkan persamaan persamaan (44)

Phase portrait of Class 1 Excitability

20 1 V (m V )

0

W nulcline

-20

-60

0

200

400

600 time (ms)

800

1000

1200

0.5 0.4

re c o v e ry v a ria ble W (m V )

0.8 -40

V nulcline

0.6

Limit cycle 0.4

W (m V )

0.2 0.3 0.2 0 0.1 0

equilibrium 0

200

400

600 time (ms)

800

1000

1200

-0.2 -50

-40

-30

-20

(a) 40

-10 0 10 membrane potential V (mV)

20

30

40

1

20

0.9

Phase portrait of Class 2 Excitability

0 0.8

W nulcline

-20 0.7 -40 -60

0

200

400

600 time (ms)

800

1000

1200

0.5 0.4

rec ov ery v ariable W (m V )

Persamaan ini digunakan untuk mencari grafik garis nol (nullclines), dan nilai akar persamaan. Selanjutnya menganalisis sistem dinamik PDB ,untuk mencari grafik garis nol dan akar-akarnya. Dengan membuat fungsi f(v,w) dan g(v,w) pada keadaan keseimbangan maka akan menjadi.

1.2

40

V (m V )

N = @( , !) = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (39) NX N = @( , !) = 0 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (40) NX

merupakan grafik garis nol umtuk dW/dt=0. Dengan melakukan simulasi menggunakan MATLAB, didapatkan grafik garis nol untuk kedua tipe propagasi 1 dan 2 dengan nilai parameter yang sama dengan simulasi sebelumnya (Iapp = tetap). Gambar 27., menampilan nulclines dengan limit cycle untuk kedua tipe.

0.6

V nulcline

Limit cycle

0.5

0.4

equilibrium

0.3 W (m V )

utama potensial membran V dan parameter pemulihan W. Analisis sistem dinamik ini meliputi pencarian titik nol (keseimbangan) dan analisis nilai dan vektor eigen untuk mengtahui karakteristik dinamik dan bifurkasi pada model. 4.2.1 Analisis linier lokal, nilai eigen dan diagram fase Dengan meninjau kembali persamaan (30) dan (31), pada keadaan keseimbangan, nilai dV/dt dan dW/dt bernilai nol. Dengan memisalkan ruas kanan pada kedua persamaan adalah f(v,w) dan g(v,w) maka persamaan (30) dan (31) menjadi.

0.3 0.2

0.2

0.1

0.1 0

0

200

400

600 time (ms)

800

1000

1200

0 -50

-40

(b)

-30

-20

-10 0 10 membrane Voltage V (mV)

20

30

Gambar 27. Diagram fase (a) tipe 1 dan (b) tipe 2 dengan Iapptetap. Untuk memahami makna kualitatif dari diagram fase tersebut, langkah selanjutnya yang dilakukan adalah mencari nilai eigen untuk menentukan jenis titik kritis (keseimbangan) pada sistem. Untuk mencari nilai eigen tersebut, maka harus dibangun suatu matrik karaktersitik yang disebut matriks jacobian (J). Dengan memasukan persamaan (30) dan (31) kedalam matriks, maka akan didapatkan, Q@(b, T) Q@(b, T) QT • ∙∙∙∙∙∙∙∙∙ (45) • = € Qb Q (b, T) Q (b, T) Qb QT

40

23 @( , !) = (− _ G _∞ ( )( − _ G ) − _w !( − _w ) − _W ( − _W ) + $_Gxx)/ ∙∙∙∙∙∙ ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (46) !∞ ( ) − ! ( , !) = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (47) '( ( )

ms/cm2, gl=2 ms/cm2gCa=4 ms/cm2 , To=1/15 s-1, VCa= 120 mV, VK=-80 mV, Vl= -60 mV, V1=-1.2 mV, V2=18 mV, V4 =17.4 mV , V3=12 mV. danIapp= 50 µA. hasil penurunan matriks dengan MATLAB didapatkan matriks (48) untuk tipe 1 dan matriks (49) untuk tipe 2 dengan nilai nilai Iapp=55 µA. dan V3=2 mV.

Dengan memasukan nilai parameter untuk tipe 1 adalah C=20 µF/cm2, gK=8 ‚=

… „ „

†

™ š –8

•‹ Œ• ’ 9ž Œ›œ •‘

ƒ( ∙∙∙∙∙∙∙∙∙∙ (48) ‚= … „ „

™ š –8

‹ Œ • ‡ˆ‰Š8 Ž 9 Œ• Œ• ’Œ“(”’./H) ‘ ‘ •‹ ••

Ÿ ¡ Œ • –ž •› •‘’ — ¢¢

.£3

•

Υ

•

/H

•‹ Œ•

Ÿ ¡ Œ • –ž •› •›’ — ¢¢

.£3

•

Υ

•

−

− cosh

‹ Œ • ‡ˆ‰Š8 Ž 9 Œ• Œ• ’Œ“(”’./H) † ‘ ‘

•‹ • ’ 9ž Œ›œ •›

ƒ( ∙∙∙∙∙∙∙∙∙∙ (49)

/H

‹ Œ Œ• Œ•

• –8 — 9

−

fungsi f(V,W) dan g(V,W) diberi masukan nilai V0 dan W0 yang dapat dicari dengan mengakarkan persamaan (41) dan (42). Pada tipe 1 dan 2, nilai akar-akar nol nya adalah, −5.375752 ªKxO 1: U H V = 8 9 !H 0.1194959

selanjutnyua pada masing-masing tipe disubstitusikan nilai V0 dan W0 pada V dan W sehingga matriks (48) dan (49) menjadi bernilai eksak. 0.4353 −29.8497 •(ªKxO 1) = 8 9 0.0009 −0.0752

−0.0331 −18.0127 • (ªKxO 2) = 8 9 0.0002 −0.1080

−

/˜ ™

−™ .

•‹ •• •

• –8 ’ 9 ™” .H •› •‘ (.£3 − /¤)( ™// ‹ Œ Œ• Œ•

• –8 — 9

− cosh

−34.968339 ªKxO 2: U H V = 8 9 !H 0.01407425

.H

.H

−

/˜ ™

−™

− ™//) − .

.

•‹ Œ• •

• –8 ’ 9 ™” ™ •› •› (.£3 − ©£)( ™//

− ™//) − .

/” ™

− 32

¥ š–8

•‹ Œ• ’ 9 Œ›œ •‘

.™

/” ™

− 32

¥ š–8

•‹ • ’ 9 Œ›œ •›

.™

¨ §∙ § ¦

¨ §∙ § ¦

Setelah didapatkan matriks jacobian, maka langkah terakhir adalah mencari nilai eigen λ pada masingmasing tipe. Hasil yang didapatkan dari simulasi MATLAB untuk nilai eigen pada tipe 1 dan 2 adalah, Y 0.3751 ªKxO 1: U . V = 8 9 Y/ −0.0149

Y −0.0705 + 0.0412K ªKxO 2: U . V = 8 9 Y/ −0.0705 − 0.0412K

dari hasil pencarian nilai eigen tersebut dapat disimpulkan bahwa titik kritis pada tipe 1 adalh titik saddle tidak stabil dengan ditandai oleh adanya nilai eigen yang betanda positif. Sedangkan nilai eigen pada tipe 2 adalah komplekskonjugat dengan suku real memiliki tanda negatif adalah titik focus yang stabil Titik kritis diatas didapatkan pada saat keadaan setimbang. Pada

24

Gambar 28. untuk diagram fase tipe 1, grafik W nulcline memotong grafik V nulclins pada 3 titik. Semua titik adalah tidak stabil. Titik kestabilan yang pertama ini merupakan tempat saat Iapp tidak cukup untuk mengeksitasi saraf sehingga akan tetap disana. Saat Iapp cukup untuk mengeksitasi, maka titik keseimbangan akan bergeser dan merubah sifat dinamiknya ke keadaan yang tidak stabil dan saraf mulai tereksitasi. 1

Biffurcation Diagram of Class 1

Gambar 29. Bifurkasi saddle-node on invariant circle (SNIC) Untuk tipe 2, memiliki jenis titik focus yang dengan diagram bifurkasi nya dapat dilihat pada Gambar 30. berikut. 1

0.8

0.9

r e c o v e r y v a r ia b le W ( m V )

0.6

Biffurcation Diagram of Class 2

0.8

Node Unstable Equilibrium 0.7

Saddle Unstable Equilibrium

Stable Focus Equilibrium re c o v e ry v a ria b le W (m V )

0.4

0.6

Threshold

Periodic Limit Cycle

0.5

0.2

0.4

no equilibrium 0

Rest State

Rest State

0.3

Periodic Limit Cycle

Excitation State

Excitation State

0.2

Excitation State -0.2 -50

-40

-30

-20

-10 0 membrane potential V (mV)

10

20

30

40

Gambar 28. Bifurkasi saddle-node pada tipe 1. pergeseran titik ini merubah jenis titik kritis node menjadi saddle Perubahan jenis titik kritis dari node menjadi saddle inilah yang merupakan suatu bifurkasi dalam sistem dinamik. Dalam hal ini nilai eigen yang bertanda positif bergerak menuju nol dan menjadi negatif sehingga menjadi stabil. Jenis bifurkasi pada tipe 1 ini adalah bifurkasi saddle-node18,19 Saat saraf memasuki keadaan eksitasi, limit cycle melewati salah satu titik kritis tidak stabil dan titik kritis tidak stabil lainnya berada di dalam nya. Sedangkan titik kritis yang stabil tidak dilewati atau berada di luar limit cycle. Jenis bifurkasi saddle-node ini adalah saddle-node on invariant circle (SNIC) bifurcation (Gambar 29).

0.1

0 -60

-50

-40

-30

-20

-10 membrane Voltage V (mV)

0

10

20

30

40

Gambar 30. Bifurkasi Andronov-Hopf pada tipe 2. Saat keadaan istirahat, tipe 2 memiliki jenis titik kritis focus stabil. Saat memasuki keadaan eksitasi, karena titik focus adalah stabil, maka ketika ada rangsangan yang cukup dari luar, saraf memulai eksitasi, jika belum cukup maka tidak akan terjadi eksitasi. Dalam hal dinamika saraf, ini berarti saraf akan mengalami eksitasi apabila ada perubahan arus terapan tertentu yang melewati nilai keadaan istirahat. Jika dilihat pada diagram bifurkasi, hanya ada 1 titik keseimbangan saja yaitu berada di dalam limit cycle.19 Oleh karena itu, daerah istirahat terletak di dalam limit cycle. Kedua sistem ini memiliki tipe bifurkasi yang berbeda. Tipe satu adalah jenis titik node yang berubah menjadi saddle saat memasuki keadaan eksitasi. Sedangkan tipe 2 adalah jenis titik focus dan tidak mengalami perubahan jenis titik kritis, namun titik kritis tersebut

25

memahami pengaruh nilai eigen dalam menjelaskan sistem dinamik pada saraf. Dengan mengganti nilai V3 pada sistem, untuk tipe 1 (Iapp=50 µA) dan 2 (Iapp=55 µA) adalah V3=18, dengan langkah yang sama pula, maka akan didapatkan nilai eigen masing-masing sebagai berikut.

kehilangan kestabilan sehingga terjadi periodic spiking. Tipe bifurkasi pada tipe 2 ini adalah bifurkasi Andronov-Hopf, seperti pada Gambar 31.

Y Y 0.3907 0.3953 U .V = 8 9 U .V = 8 9 Y/ −0.0367 −0.0346 Y/

nilai eigen pada kedua kasus adalah berlawanan tanda, sehingga kedua tipe ini memiliki jenis titik kritis saddle yang tidak stabil. Diagram bifurkasi kedua tipe dapat dilihat pada Gambar 32. Jenis titik kritis lain yang mungkin pada sistem dinamik saraf adalah titik focus. Titik ini bisa didapatkan pada kedua tipe dengan mengganti nilai V3 menjadi -3 mV. Nilai eigen masing-masing tipe akan berubah menjadi bilangan kompleks-konjugat dengan nilai masing-masing sebagai berikut,

Gambar 31. Bifurkasi Andronov-Hopf. 4.2.2 Nilai eigen dan diagram fase tipe 1 dan 2 variasi Iapp dan V3 Seperti yang telah dijelaskan sebelumnya, bahwa karakteristik sistem dinamik bergantung pada nilai inisiasi parameter yang berkaitan dengan sistem tersebut. Sebagai contoh, perubahan nilai Iapp pada persamaan akan mengubah nilai eigennya. Dengan demikian akan berubah pula karakteristik dinamiknya. Besar kecilnya perubahan parameter memiliki dua kemungkinan. Kemungkinan pertama sistem tidak akan mengubah karakteristiknya dengan jenis dan tanda nilai eigen yang tetap, namun hanya mengubah besarnya saja. Kemungkinan kedua jenis dan tanda nilai eigen akan berubah sehingga karakteristik dinamiknya akan berubah. Pada sub bab ini akan dibahas kemungkinan kedua agar lebih

Y −0.0782 + 0.0533K ªKxO 1: U . V = 8 9 Y/ −0.0782 − 0.0533K

Y −0.0720 + 0.0588K 9 ªKxO 2: U . V = 8 Y/ −0.0720 − 0.0588K

1.2

1

Calss 1 Excitability

0.9

Calss 2 Excitability

1

(a)

0.7

Initial Condition

Saddle equilibria

mem bran potenstial V (m V)

0.8

0.8

(b)

New Equilibria

0.6

Critical Point 0.4

0.6

New Equilibria Dissapear Saddle equilibria

0.5

Critical Point

0.4

0.3

0.2

0.2 0

-0.2 -50

Rest State

-40

-30

-20

-10

0

10

20

30

Initial Condition

Rest State

0.1

0 40 -50 Time (ms)

-40

-30

-20

-10

0

10

20

Gambar32. Diagram bifurkasi (a) tipe 1 dan (b) tipe 2dengan nilai V3=18 mV.

30

40

26

Maka jenis titik kedua tipe sekarang adalah nodeyang memiliki perbedaan kestabilan.Pada tipe 1 adalah tidak stabil sedangkan tipe 2 stabil. Tabel 1.Hubungan nilai V3 dan Iapp dengan bifurkasi. V3 (mV ) 12 2 18 -3

Tipe 1 (Iapp=50 mV) saddle node unstabel saddle focus

Tipe 2 (Iapp=55 mV) node stable focus saddle focus

Bifurkasi

saddle-node AndronovHopf saddle-node AndronovHopf

Dapat disimpulkan bahwa saat nilai V3 pada kedua tipe bernilai 18 mV, maka sistem tidak stabil dengan tipe bifurkasi saddle-node. Saat nilai mulai turun V3 mulai turun dan memasuki negatif (V3=-3 mV) maka sistem mulai stabil (suku real (nyata) bilangan kompleks eigen yang negatif) dan perlahan-lahan memasuki keadaan istirahat dengan tipe bifurkasi AndronovHopf. Secara menyeluruh, hubungan antara nilai parameter V3 dan Iapp dapat dilihat pada tabel 1. 4.2.3 Nilai eigen dan diagram fase tipe 1 dan 2 Iapp bergantung waktu Analisis sistem dinamik pada penjelasan sebelumnya menggunakan parameter Iapp dengan nilai yang tetap terhadap waktu. Sehingga dalam menentukan tipe bifurkasi nya agak sulit terutama dalam hal perubahan

Arus terapan DC bergantung waktu Berdasarkan persamaan (36) fungsi arus I(t) dimodelkan dengan suatu fungsi linier dengan nilai parameter α sebagai gradien laju arus terhadap waktu. Pada tipe 1 dan 2, dengan nilai α positif didapatkan bentuk propagasi seperti pada Gambar 33. Jika diperhatikan, ada tiga daerah utama pada bentuk propagasi tersebut yaitu (A) daerah pada keadaan arus mulai naik menuju keadaan eksitasi dan mulai melakukan spiking, (B) daerah saat saraf melakukan periodic spiking, dan (C) daerah berarus tinggi pada keadaan saraf tidak melakukan spiking. Karakteristik dari ketiga daerah ini berbeda dikarenakan memiliki karakteristik bifurkasi yang berbeda. 40

Class 1 Excitability

20 m e m b r a n e v o lt a g e ( m V )

Y 0.3120 Y. −0.0561 9U V = 8 9 U .V = 8 Y/ 0.0552 Y/ −0.1142

karakteristik dinamiknya. Dalam sub bab ini akan di bahas perubahan karaktersitik sistem dinamik ditinjau dari adanya perubahan nilai arus terapan terhadap waktu, apakah ada perubahan tipe bifurkasi dari keadaan istirahat ke keadaan eksitasi atau sebaliknya.nilai arus terapan bergantung waktu pada penelitian ini dibagi menjadi dua tipe berdasarkan jenis arus terapannya yaitu arus terapan DC dan AC. Pertama akan dibahas karakteristik sistem dinamik arus DC bergantung waktu, selanjutnya AC

0

(T)

-20

(A) Increasing Current State

-40

(C) Steady State

(B) Periodic Spiking State -60

0

500

1000

1500

2000

2500

time (ms) 40 20 m e m b r a n e v o lt a g e ( m V )

Pada Gambar 32, noktah merah yang memiliki label new equilibria adalah merupakan titik focus yang dimaksud. Ini dapat terjadi pada kedua tipe bahwa pada eksitasi saraf, nilai eigen akan berubah dari real (titik saddle) dan akan menghilang imaginer pada tipe biffurkasi Andronov-Hopf. Sedangkan untuk nilai Iapp dan V3 pada kedua tipe diukar yaitu untuk tipe 1 dan 2 masing-masing 50 µA, 2 mV dan 55 µA, 12 mV dan nilai eigennya adalah.

0

(A) Increasing Current State

-20

(T) (C) Steady State

-40

(B) Periodic Spiking State

-60 -80

0

200

400

600

800

1000 time (ms)

1200

1400

1600

1800

2000

Gambar 33. Tiga daerah utama propagasi (a) tipe 1 dan (b) tipe 2 dengan arus DC bergantung waktu.

27

Tabel 2. Nilai eigen masing-masing daerah pada tipe 1 dan 2 arus DC bergantung waktu. Daerah (A) t=700 ms (B) t=1500 (Transi si) t=2000 (C) t=2300

Tipe 1 0.3691 -0.0085

Titik kritis saddle

Daerah (A) t=300 ms

Tipe 2 -0.0820+0.0345i -0.0820+0.0345i

Tititk Kritis focus stable

0.3642 -0.0336 -0.0933+0.2648i -0.0933-0.2648i

saddle focus stable

(B) t=1000 (Transisi) t=1600

0.1480+0.0258i 0.1480-0.0258i -0.1079+0.2342i -0.1079-0.2342i

fucus unstable focus stable

-0.1057+0.2638i -0.1057-0.2638i

focus stable

(C) t=1800

-0.1366+0.2214i -0.1366-0.2214i

focus stable

Selanjutnya akan dibahas jenis titik kritis di tiap daerah untuk tipe 1 dan 2. Dengan nilai Imdan Iinit masingmasing pada tipe 1 dan 2 adalah 5 mV dan 10 mV dan 0 µA, didapatkan hasil dari simulasi MATLAB nilai eigen dari masing-masing daerah pada tipe 1 dan 2 yang dapat dilihat pada tabel 2. Berdasarkan data pada Tabel 2., dapat dilhat bahwa pada propagasi tipe 1 mengalami perubahan jenis titik kritis dan kestabilan. Sedangkan pada tipe 2, tidak mengalami perubahan titik kritis, hanya mengalami perubahan kestabilan saja. Pada tipe 1, dari keadaan istirahat (A) ke keadaan eksitasi (B) memiliki jenis titik kritis saddle. Seperti yang telah dijelaskan sebelumnya, meskipun tidak adanya titik node, titik ini menghilang karena sistem dalam keadaan mulai tereksitasi. Oleh karena itu, saat sistem beralih dari keadaan istirahat menuju keadaan eksitasi, jenis bifurkasi yang

terjadi adalah saddle-node. Saat memasuki daerah transisi (T), sistem mulai beralih dari keadaan eksitasi menuju istirahat. Pada tahap ini, sistem mengalami dua perubahan sekaligus yaitu perubahan jenis titik kritis dan kestabilan. Titik kritis berubah dari titik saddle tidak stabil menjadi titik focus stabil. Dari perubahan titik kritis ini dapat disimpulkan bahwa daerah transisi dari keadaan eksitasi menuju istirahat memiliki tipe bifurkasi Andronov-Hopf. Memasuiki daerah (C) yang nilai arus terapannya terlalu besar, memiliki jenis titik focus yang stabil. Jika dibadingkan dengan daerah transisi, nilai suku real memiliki nilai yang lebih besar. Ini menandakan bahwa dengan terus bertambahnya nilai Iapp, maka akan menaikan nilai eigen menuju nol dan akhirnya bertanda positif sehingga akan kembali tidak stabil. Diagram fase pada tipe 1 dapat dilihat pada Gambar 34.

28

V (m V )

50

Resting supercritical Andronov-Hopf

Saddle-Node Bifurcation

Resting

0 -50

Spiking State -100

0

500

1000

1500

time (ms)

0.5

2000

Rest State

Rest State

Stable Focus Unstable Saddle

0.4

2500

Stable Focus

Stable Node

2000

tim e (m s )

0.2

Spiking State

limit cycle

Unstable Saddle

0.3

W (m V )

2500

Stable Node

1500

1000

Spiking State

500

0.1

0 0.5

limit cycle

0

0.4

Rest State -0.1 -60

-50

-40

-30

-20

-10 V (mV)

0

10

initial condition 0.3

initial condition 20

30

W (mV)

40

40 20

Rest State

0.2

0 -20

0.1

V (mV)

-40 0

-60

Gambar 34. Diagram fase tipe 1 dengan fungsi arus DC bergantung waktu. dari keadaan istirahat menuju eksitasi.3 Sedangkan untuk daerah transisi (T), nilai eigen berubah menjadi negatif kembali sehingga sistem mulai stabil untuk memasuki keadaan istirahat. Sistem terus berosilasi dengan nilai amplitudo pulsa yang semakin melemah dan akhirnya hilang. Jenis bifurkasi yang memiliki karakteristik demikian adalah bifurkasi supercritical-Andronov-Hopf. Diagram fase untuk tipe 2 tersebut dapat dilihat pada Gambar 35.

Tipe 2 hanya memiliki satu jenis titik kritis yaitu focus. Pada propagasi ini peralihan dari (A) menuju (B) terjadi akibat perubahan kestabilan titik kritis dari stabil menjadi tidak stabil.Jika sistem tidak stabil, maka saraf akan memulai eksitasi. Jenis bifurkasi dari keadaan istirahat menuju eksitasi adalah tipe bifurkasi subcritical-Andronov-Hopf. Jenis bifurkasi ini dikatakan subcritical dikarenakan sistem megalami osilasi yang kecil saat akan melakukan transisi

V (m V )

50

0

Resting subcritical Andronov-Hopf

Resting

0 0.7

0.6

0.5

200

400

600

800

1000

Stable Focus

Rest State

1600

1800

Rest State Stable Focus

limit cycle

tim e (m s )

W (m V )

1400

2000

unstatable Focus dissapear

Spiking State 0.3

unstatable Focus dissapear

1000

limit cycle Spiking State

500

0.2

Spiking State

0.1

0 0.8

initial condition

Rest State

0.6

0

20

initial condition -0.1 -70

1200

time (ms)

1500 0.4

supercritical Andronov-Hopf

Spiking State

-50

W (mV)

0.4

0 -20 0.2

-60

-50

-40

-30

-20 V (mV)

-10

0

10

20

-40

30 0

-60

V (mV)

-80

Gambar 35. Diagram fase tipe 2 dengan fungsi arus DC bergantung waktu.

40

29

Tabel 3. Nilai eigen masing-masing daerah pada tipe 1 dan 2 arus AC bergantung waktu Daerah

Tipe 1

Max current t=150 ms Trantition t=280 Min current t=425

Titik kritis saddle

0.3770 -0.0192 0.3753 -0.0152 0.3713 -0.0105

Daerah Max current t=100ms

-0.0528+0.0487i -0.0820-0.0345i

Tititk Kritis focus stable

saddle

Trantitiont=200

saddle

Min current t=300

-0.0715+0.0407i -0.0715-0.0407i -0.0866+0.0313i -0.0866-0.0313i

focus stable focus stable

semakin menuju nilai arus minimum, kedua nilai eigen tersebut semakin mendekati angka nol. Nilai nol adalah suatu critical point yang merupakan peralihan antara keadaan stabil dan tidak stabil pada nilai eigen.20 Sedangkan untuk tipe 2, nilai eigen dari suku real bilangan kompleks-konjugat menunjukan nilai negatif yang semakin menjauhi angka nol. Ini menunjukkan bahwa sistem tersebut semakin stabil. Pada kedua tipe nilai Imax masing-masing adalah 8 µA dan 10 µA. dengan nilai ω masing-masing adalah 0.011 s-1 dan 0.016 s-1. Kedua parameter ini sangat kecil untuk mengubah karakteristik dinamik pada kedua sistem. Ini berarti nilai arus terapan AC pada model saraf adalah sangat kecil dengan tujuan untuk mengetahui bentuk propagasi saraf saja. Diagram fase untuk masing-masing tipe disajikan pada Gambar 36 dan 37.

Arus terapan AC bergantung waktu Jenis arus terapan bergantung waktu yang kedua adalah suatu arus AC yang dimodelkan sebagai suatu fungsi sinusoidal seperti pada persamaan (38). Dengan parameter ω sebagai frekuensi pulsa arus terapan. Untuk melakukan analisis sistem dinamik propagasi saraf tipe 1 dan 2, maka propagasi tersebut akan dibagi lagi menjadi beberapa daerah seperti pada analisis sebelumnya. Hasil simulasi program didapatkan jenis titik kritis pada tiap-tiap daerah disajikan dalam Tabel 3. Berdasarkan data pada Tabel 3., terlihat bahwa kedua tipe 1 dan 2 tidak mengalami perubahan jenis maupun kestabilan titik kritis. Meskipun titik tersebut tidak mengalami perubahan, tetapi sebenarnya dengan berubahnya besar nilai eigen pun akan mempengaruhi karakteristik dinamik dari sistem. Pada tipe 1 yang berjenis titik kritis saddle

Iapp (mikroA)

V (mV)

50

Tipe 2

saddle-node Bifurcation

0 -50 0 60

100

40

200

300

400

500

600

700

max current 0

100

200

300

400

500 time (ms)

600

700

Class 1

1000

800

900

1000

Class 1 1200

Iapp t=280 Iapp t=425

mobile unstable saddle

1000

mobile unstable saddle

800 time (ms)

0.6

W (mV)

900

max current

Iapp t=150

0.8

0.4

800

min current

50

Limit Cycle

400 200

Spiking State

0.2

Spiking State Limit Cycle

600

0 0.5 0

0.4

-0.2 -60

-40

-30

-20

-10 V (mV)

0

10

20 0

0.2

-20

0.1 -50

40

Rest State

0.3

Initial Condition

Rest State

20

30

40 W (mV)

-40 0

Initial Condition

-60 V (mV)

Gambar 36. Diagram fase tipe 1 dengan fungsi arus AC bergantung waktu.

30

V (mV)

50

Andronov-Hopf Bifurcation

-50 0 80 Iapp (mikroA)

Andronov-Hopf Bifurcation

0

100

200

300

40 1

400

500

700

100

800

900

1000

800

900

1000

min current

max current 0

max current 200

300

400

500 time (ms)

600

700

Class 2

0.9

Iapp t=100 Iapp t=200

0.8 0.7

Iapp t=300

1200

mobile stable focus

1000

Class 2 mobile stable focus

800 time (ms)

0.6 W (mV)

600

min current

60

Rest State 0.5

Limit Cycle

0.4

Spiking State

600

Limit Cycle 400 200

0.3

Spiking State

0 0.8

0.2

Rest State Initial Condition

0.6 0.1

Initial Condition 0 -60

-50

-40

-30

-20

-10 V (mV)

0

10

0 -20

0.2 20

30

-40 0

40 W (mV)

40 20

0.4 -60 V (mV)

Gambar 37. Diagram fase tipe 2 dengan fungsi arus AC bergantung waktu. Pada Gambar 37. dapat dilihat bahwa pada kedua tipe grafik garis nol untuk V (V nulclines) bergeser selama proses dinamik berlangsung. Keadaan grafik V nulcline yang bergerak periodik ini menyebabkan pergeseran titik keseimbangan pada sistem. Pada tipe 1 yang berjenis titik keseimbangan saddle maka akan bergerak naik turun mengikuti V nulcline yang berosilasi. Begitu pula untuk tipe 2 yang berosilasi pula.Catatan bahwa titik kritis pada kedua sistem tidak mengalami perubahan jenis maupun kestabilan selama berosilasi. Jika nilai Imax dan ω divariasikan dengan interval nilai yang cukup besar, akan ada dua kemungkinan bahwa sistem akan mengubah jenis dan kestabilan titik kritis karena nilai arus terapan dengan fluktuasi yang tinggi, atau sistem tetap mempertahankan karakteristik dinamik awal nya (tidak mengalami perubahan karakteristik titik kritis). Kedua kemungkinan ini tidak dibahas pada penelitian ini karena dalam analisis arus terapan AC ini sudah cukup untuk mengetahui karakteristik dinamik suatu propagasi dengan arus terapan yang sangat kecil. 4.3 Solusi Numerik Propagasi Saraf Terkopel Model jaringan saraf yang dibahas sebelumnya merupakan hasil model jaringan saraf yang diwakili oleh satu sel tunggal. Jaringan saraf merupakan suatu gabungan fungsional

dari banyak saraf dengan sifat dan karakteristik tertentu. Dengan demikian dalam penelitian ini dibangun suatu model saraf kompleks yang melibatkan banyak saraf yang saling terhubung secara fungsional. Solusi dari model yang dibangun menganggap bahwa saraf terhubung satu dengan yang lainnya secara sinaptik. Kata sinaptik ini berasal dari salah satu komponen sel saraf pada ujung bagian akson yang terhubung dengan badan sel lainnya disebut synapses . Melalui bagian inilah sel satu dengan yang lainnya bertukar informasi.12 Hubungan sinaptik ini memiliki sifat tertentu dilihat dari bagaimana hubungan tersebut terjadi pada dua sel saraf yang terkopel.17 • Sinaptik elektrik: merupakan suatu pengiriman informasi dari satu sel ke sel lain berdasarkan peristiwa difusi linier pada potensial membran saraf terkopel. • Sinaptik kimia: merupakan suatu pengiriman informasi secara nonlinier yang melibatkan fenomena sinkronisasi pada model saraf pemacu (excitatory) dan penghambat (inhibitory).21 Pada penelitian ini dibahas tipe sinaptik elektrik. Model saraf yang dibangun pada penelitian ini adalah suatu model saraf dengan asumsi bahwa suatu jaringan saraf kompleks dapat dimodelkan oleh dua saraf terkopel yang saling terhubung

31

secara sinaptik.19,21 Jika bahasan mengenai dua saraf terkopel ini dapat dijelaskan, maka akan mudah membangun sistem banyak saraf yang saling terkopel satu dengan lainnya secara sinaptik. 4.3.1 Model saraf terkopel Model yang digunakan pada penelitian ini merupakan suatu model saraf terkopel hasil penggabungan dan modifkasi dari model saraf terkopel sebelumnya, sehingga model yang dipakai pada simulasi merupakan suatu model saraf sinaptik terkopel. Secara umum, model untuk banyak saraf telah dipublikasikan oleh Hoppensteadt dan Izhikevich (1997) dengan hanya memperhatikan kopling potensial membran antar sel3,4 seperti pada persamaan (28) dan (29). Jika persamaan (29) digabungkan dengan persamaan (28), maka akan menjadi. ′

=

@(b ) + q ∑ h. -

.

.—®¯°8σj±² ’³² m9

´ (50)

Model pada persamaan (50) merupakan suatu model dengan mengasumsikan bahwa semua sel saraf dalam suatu sistem adalah saling terkopel dan tidak memperhatikan nilai potensial pembalik setelah melakukan kopling dengan sel saraf lain. Oleh karena itu diusulkan suatu model yang menambahkan pengaturan nilai potensial pembalik dan keterhubungan antar sel.17 = @(b ) − (b −

µ )q g

fℎ g h.

g jbg m

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (51) K = 1,2,3, … . . , os = 1,2,3, … … o

ℎg

1, sK·G K NGo s XOMℎSISo =¶ ¸∙ 0, sK·G K NGo s XG· XOMℎSISo ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (52)

Jika persamaan (51) disederhanakan dengan asumsi bahwa kopling antar sel saraf dipengaruhi oleh suatu arus sinaptik, maka fungsi sinaptik kopling dapat dibentuk sebagai fungsi potensial membran tiap saraf ditambah dengan fungsi arus sinaptik Isyn.. ′

g = @(b ) + $µ¹ ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (52)

$µ¹ = −(b − g

µ )q f ℎ g h.

g jbg m

K = s = 1,2,3, … . . , o

∙∙ (53)

Sekarang persamaan (52) dan (53) akan ditransformasi ke dalam model ML. Dengan mensubstitusikan persamaan (34) dan (35) kedalam fungsi f(vi), makadidapatkan persamaan berikut.

K = 1,2,3, … . . , os = 1,2,3, … … o

′

antara kedua saraf terhubung atau tidak, dengan ketentuan sebagai berikut.

∙

Vs merupakan potensial pembalik dengan anggapan bahwa pada hubungan sinaptik kimia, hubungan sinaptik ini merupakan jenis penghambat (inhibitory). Sedangkan hij merupakan suatu parameter Heaviside yang menentukan apakah

N =− NX

∞(

− −

)(

−

)

!( − ) g " ( − " ) + $ %%

+ $µ¹ ( g ) ∙∙∙∙∙∙∙∙∙∙∙∙ (54) N! !∞ ( ) − ! = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (55) NX '( ( ) K = s = 1,2,3, … . . , o g

Persamaan inilah yang merupakan model sinaptik kopling Morris-Lecar dengan nilai arus terapan yang dapat divariasikan. Untuk model kopling 2 saraf dengan nilai n=2, maka model kopling menjadi. N . . ) =− ∞ ( . )( . − NX − !. ( . − ) − " ( . − " ) + $./%% ./ + $µ¹ ( / ) ∙∙∙∙∙∙∙∙ (56. G) .( ) N!. !∞ . − !. = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (56. I) .( ) NX '( .

32 N / =− NX

/ ∞ ( / )( /

)

−

Hasil simulasi untuk tipe 1didapatkan hasil seperti pada Gambar 38.

− !/ ( / − ) − " ( / − " ) + $ /.%% /. + $µ¹ ( . ) ∙∙∙∙∙∙∙∙∙ (56. P) /( ) N!/ !∞ / − !/ = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (56. N) /( ) NX '( /

40

neuron 1

V (mV)

20

0

-20

-40

-60

0 40

V (mV)

100

200

300

400

500

600

700

800

900

1000

200

300

400

500 time (ms)

600

700

800

900

1000

neuron 2

20

0

-20

-40

-60

0

100

Gambar 38. Dua saraf tipe 1 nonkopling.εi,2=0.

Parameter Iapp dapat berupa arus DC tetapmaupun bergantung waktu atau arus AC. Dalam penelitian ini, akan dibahas jenis Iapp sebagai konstanta dan bergantung waktu AC. Untuk DC tidak akan dibahas. 4.3.2 Solusi numerik model saraf terkopel Iapp tetap. Agar memahami lebih lanjut fenomena kopling ini, dengan meninjau kembali persamaan (8), parameter ø merupakan skala waktu propagasi. Secara garis besar, parameter inilah yang menyebabkan perbedaan fase pada propagasi saraf. Dengan memisalkan dua buah sel saraf dengan tipe eksitasi yang sama yaitu keduanya merupakan tipe 1, atau keduanya merupakan tipe 2. Maka akan dibuat simulasi sinkronisasi kedua saraf tersebut dengan nilai ø yang sama, atau berbeda. Pada kasus pertama dengan nilai ø yang sama yaitu ø=1/15 s-1 , pada keadaan terisolasi (bebas tidak saling mempengaruhi) εi,2=0, hij=0, dan nilai parameter Vs= 2 mV, σ=0.01, θ=-40 mV.

Berdasarkan hasil yang didapat pada Gambar 38., dapat dilihat bahwa dengan nilai ø yang sama, skala waktu propagasi kedua saraf sama. Yang membedakan hanya fase awal nya saja, pada saraf 1 memiliki nilai potensial awal adalah -40 mV, sedangkan saraf kedua 0 mV. Kemudian kedua saraf dikopelkan (εi,2≠0, hij=1), dengan mengubah εi= 0.5 mS/cm2 dan ε2=1.25 mS/cm2, maka kedua saraf sudah terkopel, dan didapatkan hasil seperti pada Gambar 39. Dapat dilihat bahwa baik fase maupun frekuensi eksitasi sudah berbeda. Ini terjadi dikarenakan propagasi masing-masing saraf dipengaruhi satu sama lain dengan kekuatan kopling yang berbeda (εi ≠ ε2) sehingga menghasilkan propagasi yang berbeda. Untuk mensinkronkan propagasi kedua saraf tersebut, maka kekuatan kopling antara kedua saraf tersebut harus

V (mV)

50 0

V (mV)

-50 0 50

neuron 1 100

200

300

400

500

600

700

800

900

1000

200

300

400

500 time (ms)

600

700

800

900

1000

0 -50

neuron 2 0

100

0.45 0.4

Class 1 Excitability

1200

0.35

1000 800 time (ms)

W (mV)

0.3 0.25 0.2

600 400 200

0.15

0 0.5

0.1

0.4

40 0.3

0.05

20 0

0.2 0 -50

-20

0.1 -40

-30

-20

-10 0 V (mV)

10

20

30

40 W (mV)

-40 0

-60 V (mV)

Gambar 39. Tipe 1 dua saraf terkopel (εi,2≠ 0, hij=1), non-sinkronisasi.

33

50 V (m V)

0 -50

V (m V)

-100 50 0

neuron 1 100

200

300

400

500

600

700

800

900

1000

200

300

400

500 time (ms)

600

700

800

900

1000

0 -50

neuron 2 0

100

0.5

Class 1 Excitability

0.45

1200 0.4 1000 0.35 tim e (m s)

800 W (m V )

0.3 0.25

600 400

0.2 200 0.15

0 0.5

0.1

0.4

40 0.3

0.05

20 0

0.2 0 -60

-20

0.1 -50

-40

-30

-20

-10 V (mV)

0

10

20

30

40

0

W (mV)

-40 -60 V (mV)

Gambar 40. Tipe 1 dua saraf terkopel (εi,2≠ 0, hij=1), tersinkronisasi (εi = ε2=0.5 mS/cm2). diseragamkan (εi = ε2). Dengan mengubah nilai εi = ε2 = 0.5 mS/cm2 dan Vi=V2=0, maka didapatkan propagasi saraf tersinkronisasi seperti pada Gambar 40. Meskipun propagasi yang terjadi memiliki fase yang berbeda, namun kedua saraf memiliki skala waktu perambatan yang sama seperti dilihat pada ruang fase pada Gambar 40. kedua saraf yang saling berhimpitan.

0.5

non-coupled, eps1=eps2=0

Class 1 Excitability

0.4

W (mV)

Kasus kedua pada keadaan kedua saraf memiliki nilai skala waktu ø yang berbeda. Dengan keadaan yang sama seperti pada keadaan sebelumnya sedangkan nilai ø untuk masing-masing saraf adalah ø1=1/15 s-1dan ø2= 1/20 s-1, didapatkan diagram fase saraf seperti pada gambar 41. Dapat disimpulkan bahwa pada propagasi tipe 1 model dua saraf terkopel dengan skala waktu yang berbeda, sinkronisasi sangat sulit dilakukan. Hasil 0.5

No Synchronization

0.4

0.3

0.3

Near-Synchronization

0.2

0.2

0.1

0.1

Far-Synchronization 0 -60

-40

-20

0

20

0 -50

40

V (mV)

W (mV)

0.5

coupled, eps1=0.25 eps2=1.25

0.5

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

-20

0

20

40

0 -60

-40

-20

0

20

-30

-20

-10

0

10

20

30

40

0.5

0.4

-40

-40

coupled, eps1=0.85 eps2=1.25

coupled, eps1=0.5 eps2=0.5

0.4

0 -60

coupled, eps1=1.5 eps2=2.0

40

0 -60

-40

-20

0

20

V (mV)

Gambar 41. Sinkronisasi kopling saraf tipe 1 dengan nilai skala waktu berbeda.

40

34

V (mV)

50 0 -50 V (mV)

50

neuron 1 0

100

200

300

400

500

600

700

800

900

1000

200

300

400

500 time (ms)

600

700

800

900

1000

0 -50

neuron 2 0

100

0.5 0.45

1200 0.4 1000 0.35 time (ms)

800 W (mV)

0.3 0.25

600 400

0.2 200 0.15 0 0.5

0.1

0.4

0.05 0 -50

40 0.3

Class 2 Excitability

20 0

0.2

-20

0.1 -40

-30

-20

-10 0 V (mV)

10

20

30

40

-40 0

W (mV)

-60 V (mV)

Gambar 42. Tipe 2 dua saraf terkopel (εi,2 ≠ 0, hij=1), non-sinkronisasi. Nilai parameter yang dipakai adalah ε1=0.25 mS/cm2, dan ε2=1.25 mS/cm2 ,Iapp= 55 µA, V3= 2 mV. Sedangkan untuk mensinkronkan dua saraf terkopel tersebut, maka ε pada kedua saraf diseragamkan menjadi ε1= ε2=0.5 mS/cm2. Didapatkan hasil seperti pada Gambar 43. Untuk sinkronisasi dengan skala waktu yang berbeda (ø1=1/15 s-1dan ø2= 1/20 s-1) baik tipe 1 dan 2 sangat sulit dilakukan. Pada tipe 2, untuk mendekati sinkronisasi, nilai εi dan ε2 masingmasing adalah 2 mS/cm2 dan 2.5 mS/cm2. Hasil variasi nilai ε lainnya dapat dilihat pada Gambar 44.

yang didapatkan hanya mendekati sinkronisasi tapi belum tersinkronisasi. Harus diperhatikan bahwa parameter skala waktu adalah tidak sama dengan beda fase propagasi antara kedua saraf. Jika dua saraf memiliki perbedaan fase propagasi, maka akan lebih mudah tersinkronisasi dibandingkan dengan dua saraf yang berbeda skala waktu propagasinya. Kesimpulan tersebut diambil berdasarkan hasil yang terlihat pada Gambar 40 dan 41. Pada propagasi tipe 2, didapatkan hasil kopling saraf seperti pada Gambar 42.

50 V (mV)

0 -50

V (m V)

-100 0 50

neuron 1 100

200

300

400

500

600

700

800

900

1000

200

300

400

500 time (ms)

600

700

800

900

1000

0 -50

neuron 2 0

100

0.5 0.45

1200 0.4 1000 0.35 time (m s)

800 W (m V)

0.3 0.25

600 400

0.2 200 0.15

0 0.5

0.1

0.4

0.05 0 -60

40 0.3

Class 2 Excitability

20 0

0.2

-20

0.1 -50

-40

-30

-20

-10 V (mV)

0

10

20

30

40 W (mV)

0

-40 -60 V (mV)

Gambar 43. Tipe 2 Dua saraf terkopel (εi,2 ≠ 0, hij=1), tersinkronisasi (εi = ε2=0.5 mS/cm2).

35

0.5

non-coupled, eps1=eps2=0

Class 2 Excitability

W (m V)

0.4

0.5

No Synchronization

0.4

0.3

0.3

Near-Synchronization

0.2

0.2

0.1

0.1

Far-Synchronization 0 -60

-40

-20

0

20

0 -50

40

W (mV)

V (mV) 0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

-40

-20

0

20

40

0 -60

-40

-20

0

20

-40

-30

-20

-10

0

10

20

30

40

coupled, eps1=1.75 eps2=2.5

coupled, eps1=0.5 eps2=0.5

coupled, eps1=0.25 eps2=1.25 0.5

0 -60

coupled, eps1=2 eps2=2.5

0.5

40

0 -60

-40

-20

0

20

40

V (mV)

Gambar 44. Sinkronisasi kopling saraf tipe 2 dengan nilai skala waktu berbeda. Propagasi pada tipe 1 dan 2 diatas hanya melibatkan nilai arus Iapp tetap dan kekuatan kopel antar kedua saraf. Sedangkan untuk perbedaan skala waktu menyebabkan dua saraf terkopel sangat sulit untuk tersinkronisasi. Jika nilai parameter lain ikut divariasikan seperti potensial pembalik Vs ,jenis kopling menjadi suatu saraf pemacu excitatory, dan nilai laju kopling σ, maka akan didapatkan hasil yang lebih bervariasi dari hasil simulasi diatas. Dengan demikian fenomena sinkronisasi ini sangat bergantung dengan karakteristik propagasi tiap-tiap saraf dalam suatu jaringan kompleks. 4.3.3 Solusi numerik model saraf terkopel Iapp AC bergantung waktu. Seperti telah yang dijelaskan sebelumnya mengenai bahasan pengaruh arus Iapp bergantung waktu yang akan dibahas adalah merupakan fungsi arus AC. Dengan mensubstitusikan persamaan (38) ke dalam persamaan (54) dan (55). N =− NX

∞(

)(

−

)

− !( − ) − "( − ") g g + $z { sin (~X) + $ a g + $µ¹ ( g ) ∙∙∙∙∙∙∙∙∙∙∙∙ (57)

N! !∞ ( ) − ! = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (58) NX '( ( ) K = s = 1,2,3, … . . , o

Pada sistem kopling dua saraf, maka persamaan diatas menjadi. N . =− NX

. ∞ ( . )( .

−

)

− !. ( . − ) − "( . − ") . . + $z { sin (~X) + $ a ./ + $µ¹ ( / ) ∙∙∙∙∙∙∙∙ (59. G) .( ) N!. !∞ . − !. = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (59. I) .( ) NX '( . N / / ) =− ∞ ( / )( / − NX − !/ ( / − ) − "( / − ") / / + $z { sin (~X) + $ a /. + $µ¹ ( . ) ∙∙∙∙∙∙∙∙∙ (59. P) /( ) N!/ !∞ / − !/ = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (59. N) /( ) NX '( /

Hasil simulasi untuk tipe 1 dengan variasi nilai ε dihasilkan pada Gambar 45

36

Class 1 Excitability

Injected AC current 80 60 40

0

100

200

300

400

500

600

700

800

900

1000

200

300

400

500

600

700

800

900

1000

400

500

600

700

800

900

1000

500

600

700

800

900

1000

500 time (ms)

600

700

800

900

1000

Not coupled hij=0

100 0 -100

0

100

coupled hij=1, eps1=0.25, eps2=1.25, not synchron 100 0 -100

0

100

200

300

coupled hij=1, eps1=0.5, eps2=0.5, synchron with different phase

50 0 -50

0

100

200

300

400

coupled hij=1, eps1=0.5, eps2=0.5, synchron with same phase 50

0

-50

0

100

200

300

400

Gambar 45. Sinkronisasi tipe 1 dengan arus terapan AC. dilihat dari hasil simulai bahwa kedua saraf memiliki propagasi yang identik dengan beda fase -40 mV. Untuk membuat fase kedua saraf sama, maka nilai awal potensial membran kedua saraf diseragamkan menjadi 0 mV. Didapatkan propagasi yang identik. Ini berarti kedua saraf telah sinkron dengan fase propagasi yang sama. Dengan mengganti nilai Imax= 10 mV, Iinit=55 mA, dan V3=2 mV, untuk kedua saraf pada tipe 2, maka didapatkan hasil seperti pada Gambar 46. Hasil yang unik didapatkan pada simulasi tipe 2. Pada keadaan saraf tidak terkopel, kedua saraf diberi fase yang berbeda -40 mV, propagasi kedua saraf hanya berbeda pada bagian awal yaitu saat t<400 ms. Saat t>400 ms, fase kedua saraf hampir sama dan akhirnya sefase. Terlihat pada grafik garis hijau dan merah yang saling berhimpitan. Keadaan ini dapat terjadi walaupun kondisi kedua saraf tidak terkopel (terisolasi satu sama lain). Sedangkan saat keadaan saraf

Hasil simulasi pada tipe 1 menunjukan bahwa saat nilai arus menuju positif, maka frekuensi eksitasi lebih tinggi dibandingkan dengan saat nilai arus pada keadaan negatif. Kondisi ini berlaku pada keadaan terkopel maupun tidak. Hal yang membedakan adalah saat kondisi tidak terkopel, kedua saraf memiliki frekuensi yang lebih lebar dibandingkan saat keadaan terkopel. Sama seperti sebelumya, hal ini disebabkan oleh adanya potensial batasan yang memiliki tipe inhibitory. Ini berarti bahwa saat potensial saraf pertama mencapai nilai potensial pembalik Vs=40 mV, maka seolah saraf dua akan kembali menaikan nilai potensial saraf tersebut, sehingga akan lebih cepat tereksitasi, begitu juga sebaliknya.18 Saat keadaan saraf mulai terkopel, dapat dilihat bahwa ketika nilai ε1≠ε2 maka saraf belum tersinkronisasi. Saat nilai ε1=ε2 ,maka kedua saraf sudah sinkron, namun masih memiliki perbedaan fase propagasi. Ini dapat

Class 2 Excitability

Injected AC current 80 60 40

0

100

200

300

400

500

600

700

800

900

1000

200

300

400

500

600

700

800

900

1000

400

500

600

700

800

900

1000

500

600

700

800

900

1000

500 time (ms)

600

700

800

900

1000

Not coupled hij=0

100 0 -100

0

100

coupled hij=1, eps1=0.25, eps2=1.25, not synchron 50 0 -50 50

0

100

200

300

coupled hij=1, eps1=0.5, eps2=0.5, synchron with different phase

0 -50

0

100

200

300

400

coupled hij=1, eps1=0.5, eps2=0.5, synchron with same phase 50

0

-50

0

100

200

300

400

Gambar 46. Sinkronisasi tipe 2 dengan arus terapan AC.

37

mulai terkopel dengan kekuatan yang berbeda, propagasi kedua saraf tidak sama. Saat nilai ε bernilai sama dan sefase, maka kedua saraf telah tersinkronisasi. Keadaan ini dinamakan keadaan sinkronisasi fase terkunci phase locking Synchronization.22 Perbedaan antara propagasi tipe 1 dan 2 saraf terkopel hampir sama dengan pada saat tidak terkopel. Untuk saraf terkopel tipe 1, saat nilai arus memasuki negatif, maka frekuensi spike akan menurun. Sedangkan untuk tipe 2, saat nilai arus negatif, maka tidak akan terjadi spike, melainkan terjadi pemuluran waktu delay yang menyebabkan bursting. Khusus untuk tipe 2, antara keadaan terkopel dan tidak adalah saat terkopel, frekuensi bursting akan lebih cepat terjadi dibandingkan saat tidak terkopel. Ini berkaitan dengan penjelasan sebelumnya pada tipe 1, bahwa jenis kopling diatas adalah merupakan jenis inhibitory. 4.4 Solusi Numerik pada n Saraf Terkopel Agar lebih memahami konsep mengenai model kopel saraf, dan untuk mendekati kenyataan sesungguhnya bahwa jaringan saraf merupakan suatu sistem yang kompleks, maka khusus pada sub bab ini akan ditambahkan bahasan mengenai sistem kopling saraf dengan jumlah lebih dari 2 sel saraf. Telah dijelaskan sebelumnya bahwa model 2 saraf yang saling terkopel merupakan representasi dari suatu jaringan kompleks pada saraf. Model kopel 2 saraf ini merupakan dasar pemikiran bahwa sistem saraf kompleks merupakan susunan atas banyak sistem dua saraf terkopel yang saling berhubungan. Dengan demikian, penjelasan mengenai n saraf pada sistem saraf akan dapat dijelaskan dengan sistem dua saraf terkopel. Pembahasan pada sub bab ini yaitu untuk sistem saraf dengan kopel n=2,3, dan 4 dengan arus terapan AC. Untuk n=2, telah dibahas sebelumnya, sedangkan untuk n>4, tidak akan dibahas dengan asumsi bahwa bahasan mengenai n=2,3, dan 4 sudah

dapat mewakili fenomena sinkronisasi pada sistem saraf terkopel. 4.4.1 Solusi numerik pada 3 saraf terkopel Model umum kopling n buah saraf seperti pada persamaan (54) dan (55), pada kopel n=3, maka model kopling tiga saraf akan menjadi. N . =− NX

. ∞ ( . )( .

−

)

− !. ( . − ) − "( . − ") . . + $z { sin (~X) + $ a . + $µ¹ ( / , 2 ) ∙∙∙∙ (60. G) .( ) N!. !∞ . − !. = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (60. I) .( ) NX '( . N / / ) =− ∞ ( / )( / − NX − !/ ( / − ) − "( / − ") / / + $z { sin(~X) + $ a / ( + $µ¹ . , 2 ) ∙∙∙∙ (60. P) /( ) N!/ !∞ / − !/ = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (60. N) /( ) NX '( / N 2 2 ) =− ∞ ( 2 )( 2 − NX − !2 ( 2 − ) − "( 2 − ") 2 2 + $z { sin(~X) + $ a 2 ( + $µ¹ . , / ) ∙∙∙∙ (60. O) 2( ) N!2 !∞ 2 − !2 = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (60. @) 2( ) NX '( 2

Dengan fungsi Isyn adalah. . $µ¹ = −(b.

−

µ ) -ℎ./ q./ º

1

» 1 + expjσ(b/ − t. )m 1 + ℎ.2 q.2 º »´ (61. G) 1 + expjσ(b2 − t. )m / $µ¹ = −(b/ 1 − µ ) -ℎ/. q/. º » 1 + expjσ(b. − t/ )m 1 »´ (61. I) + ℎ/2 q/2 º 1 + expjσ(b2 − t/ )m

38 2 $µ¹ = −(b2

−

Sehingga didapatkan model kopling 3 saraf yang dapat divariasikan kekuatan kopling dan keterhubungannya dengan parameter εij dengan hij. Untuk memahami konektifitas antar saraf pada sistem kopling tiga saraf ini, maka kan diilustrasikan sutau diagram mengenai kopling tersebut antara saraf 1 dan 2, saling terkopel dan sinkron. Sedangkan antara saraf 2 dan 3 terkopel, namun tidak sinkron. Beda hal dengan hubungan saraf 1 dan 3 yaitu dengan nilai h13=0, yang ditandai dengan tidak adanya garis penghubung dari saraf 1 ke 3, berarti ini tidak terhubung. Sedangkan untuk h31=1 dengan garis putus-putus dari saraf 3 ke 1, menunjukan bahwa saraf 3 terhubung dengan saraf 1 namun tidak sinkron. Untuk mencapai suatu keadaan sinkronisasi, maka antara kedua saraf harus saling terkopel.23 Hasil simulasi untuk keadaan Gambar 47., dengan nilai ε12= ε21=0.5, ε23=0.25, ε32=1.25, ε13=0, ε31=0.5 mS/cm 2 untuk Gambar 47.b, dan εij=0.5 untuk Gambar 47.a, dengan nilai arus terapan AC yang sama, maka didapatkan hasil simulasi untuk propagasi tipe 1 seperti pada gambar 48.

1

µ ) -ℎ2. q2. º

» 1 + expjσ(b. − t2 )m 1 + ℎ2/ q2/ º »´ (61. P) 1 + expjσ(b/ − t2 )m seperti pada Gambar 47. coupled n synchron coupled not synchron uncoupled

neuron 1

(a) neuron 2

neuron 3

coupled n synchron coupled not synchron uncoupled

neuron 1

(b) neuron 2

neuron 3

Gambar 47. Sistem 3 saraf terkopel.(a) terkopel dan sinkron(b) kopel tidak sempurna. Ilustrasi pada Gambar 27 diatas menjelaskan bagaimana salah satu kemungkinan keadaan ketika ketiga saraf tersebut terhubung. Dapat dilihat bahwa 3 Coupled Class 1 Excitability 80 60 40

neuron 1

neuron 3

neuron 2

Injected AC current 0

100

Not coupled hij=0

(a) 200

300

400

500

600

700

800

900

1000

100

(b)

0 -100

0

100

200

300

400

500

600

700

800

900

1000

400

500

600

700

800

900

1000

coupled hij=1, epsij=0.5, not same phase

100

(c)

0 -100

0 full 50

100

200

coupled hij=1, epsij=0.5, synchron

300

(d)

0 -50

0

100

not perfect coupled

200

300

400

500

600

700

800

900

1000

100

(e)

0 -100

0

100

200

300

400

500 time (ms)

600

700

800

900

1000

Gambar 48. Sinkronisasi 3 saraf terkopel tipe 1 dengan variasi kemungkinan keadaan terkopel. (a) arus terapan AC. (b) tidak terkopel. (c) terkopel dengan fase berbeda. (d) terkopel dan tersinkronisasi. (e) kopel tidak sempurna

39

Propagasi tiga saraf terkopel menunjukan bahwa pada keadaan saling terisolasi, saraf 1, 2, maupun tiga saling tidak mempengaruhi. Sedangkan pada Gambar 28.b menunjukan bahwa antara saraf 2 dan 3 dapat tersinkronisasi walaupun bukan dalam keadaan kopel. Ini telah dibahas pada sub bab sebelumnya terakit dengan karakteristik dinamik tiap-tiap saraf. Sedangkan ketika ketiga saraf saling terkopel, namun tidak tersinkronisasi, maka ketiga saraf ini menempati kedudukan different phase synchronization (dps) (Gambar 48.(c)), yang berarti sistem sudah tersinkronisasi namun memiliki fase propagasi yang berbeda.22 Untuk Gambar 48.d, sistem tersinkronisasi sempurna dengan fase propagasi yang sama. Keadaan ini disebut same phase synchronization (sps) yang berarti bahwa sistem tersinkronisasi dengan fase propagasi yang sama. Sedangkan untuk Gambar 48.e merupakan hasil propagasi saraf dari ilustrasi pada Gambar 47.b. dapat dilihat bahwa antara saraf 1 (garis merah) dan 2 (garis hijau) tersinkronisasi dengan fase yang berbeda. Pada saraf 2 dan 3 (garis biru), terkopel namun tidak tersinkronisasi karena ε23≠ε32. Antara saraf 1 dan 3 tidak memiliki hubungan kopel sempurna dan sama sekali tak tersinkronisasi. Hasil yang didapatkan ini, mengasumsikan bahwa ketiga saraf pada saat awal propagasinya memiliki fase yang berbeda, yaitu masing-masing nilai awal potensial pada saraf 1, 2, dan 3 adalah -60 mV, -30 m, dan 0 mV. Berbagai variasi kondisi kopling pada sistem 3 saraf dapat dilakukan dengan terlebih dahulu membuat suatu diagram sistem kopling seperti pada Gambar 47. Dengan diagram tersebut, dapat menjelaskan apakah tiap-tiap saraf terhubung (hij≠0), apakah sinkron (εi=εj), dan bagaimana kekuatan kopling tersebut. Untuk lebih memahaminya, diberikan dua contoh diagram 3 kopling tersebut dengan keadaan yang berbeda.

neuron 1

coupled n synchron coupled not synchron uncoupled

(a)

neuron 2

neuron 3

neuron 1

coupled n synchron coupled not synchron uncoupled

(b)

neuron 2

neuron 3

Gambar 49. Sistem 3 saraf terkopel. (a)variasi 1. (b) variasi 2. Berdasarkan Gambar 49., pada Gambar 49.a menunjukan kemungkinan ketika saraf 1 dan 2 tidak terhubung sama sekali. (h12=0), saraf 2 dan 3 terkopel (h23≠0) dan tersinkronisasi (ε2=ε3). Antara saraf 1 dan 3 tidak terkopel sempurna (h31=1,h13=0). Sedangkan untuk Gambar (b) merupakan sistem dengan kopel tidak sempurna (hij≠hji). Sebagai contoh, h12=1, h21=0. Hasil sistem diagram pada Gambar 49.dengan nilai εij untuk semua kopling adalah 0.5diberikan pada Gambar 50. 3 Coupled Class 1 Excitability

neuron 1

neuron 3

neuron 2

variation (a) different phase

100

(a)

0

-100

0

100

200

300

400

500

600

700

800

900

1000

variation (a) same phase

100

(b)

0

-100

0

100

200

300

400

500

600

700

800

900

(c)

0

-100

1000

variation (b) different phase

100

0

50

100

variation (b) same phase

200

300

400

500

600

700

800

900

(d)

0

-50

1000

0

100

200

300

400

500 time (ms)

600

700

800

900

1000

Gambar 50. Sinkronisasi 3 saraf terkopel tipe 1 dengan variasi kemungkinan keadaan terkopel. (a) variasi 1 beda fase (b) variasi 1 sefase (c) varaisi 2 beda fase (d) varaisi 2 sefase. Berdasarkan Gambar 50. dapat dilihat bahwa untuk variasi 1 hanya saraf

40

2 dan 3 yang tersinkronisasi walaupun belum sempurna. Ini ditandai dengan garis hijau dan biru yang berdekatan dengan frekuensi propagasi yang sama. Sedangkan untuk variasi 2, masingmasing saraf saling terkopel namun tidak sempurna. Berdasarkan hasil yang didapat dengan nilai kekuatan kopel εijadalah sama untuk setiap saraf adalah 0.5 mS/cm2, maka tetap dapat terjadi sinkronisasi pada sistem tersebut dengan ditandai oleh propagasi yang saling berhimpitan (Gambar 50.d). Dapat disimpulkan bahwa untuk mencapai keadaan phase locking baik pada keadaan dps maupun sps, terutama akan ditentukan oleh kekuatan kopling yang diwakili parameter εij dibandingkan dengan konektifitas nya apakah saling terhubung (hij=1)atau tidak (hij=0). Pada kopling propagasi tipe 2, dengan diagram yang sama seperti pada Gambar 47 didapatkan hasil seperti pada gambar 51. 3 Coupled Class 2 Excitability 80

neuron 3

neuron 2

0

100

Not coupled hij=0

200

300

400

500

600

700

800

900

40

0

-60

200

300

400

500

600

700

800

900

1000

0

(b)

-20

0

100

200

300

400

500

600

700

800

900

1000 -40 -60

0

100

200

full coupled hij=1, epsij=0.5, synchron

300

400

500

600

700

800

900

1000

(d)

0

0

100 not perfect coupl ed

200

300

400

500

600

700

800

900

(e)

0

0

100

200

300

400

500 time (ms)

600

700

800

900

0

100

200

300

400

500 time (ms)

600

700

800

900

1000

Gambar 52. Sinkronisasi 3 saraf terkopel tipe 2 dengan variasi kemungkinan keadaan terkopel (a) variasi 2 beda fase (b) variasi 2 sefase.

1000

100

-100

100

20

(b)

(c)

-50

0

variation (b) same phase

coupled hij=1, epsij=0.5, not same phase

50

(a)

20

40

0 -100

neuron 3

neuron 2

60

1000

0

100

neuron 1

variation (b) different phase

-40

(a)

100

-100

3 Coupled Class 2 Excitability

-20

Injected AC current

60 40

neuron 1

mulai tersinkroisasi dengan kondisi dps (c). Saat fase ketiga saraf adalah sama, maka tercapai keadaan sps (d).sedangkan untuk Gambar (e), hanya saraf 1 (garis merah) dan 2 (garis hijau) yang tersinkronisasi. Sinkronisasi yang terjadi antara keduanya adalah dps. Sedangkan untuk saraf 3 tidak tersinkronisasi baik dengan saraf 1 maupun 2, namun terkopel dengan saraf 2 dengan kekuatan kopling yang berbeda (ε23=0.25, ε32=1.25). Untuk membuktikan pernyataan sebelumnya yang menjelaskan bahwa sinkronisasi lebih ditentukan oleh kekuatan kopling εij dibandingkan dengan konektifitas nya hij, maka diagram pada Gambar 49.b akan disimulasikan untuk propagasi tipe 2. Hasil yang didapatkan disajikan pada gambar 52.

1000

Gambar 51. Sinkronisasi 3 saraf terkopel tipe 2 dengan variasi kemungkinan keadaan terkopel. (a) arus terapan AC (b) tidak terkopel (c) terkopel dengan fase berbeda (d) terkopel dan tersinkronisasi (e) kopel tidak sempurna. Saat ketiga saraf tidak terkopel, maka fenomena bursting masih terlihat. Saat ketiga saraf terkopel dengan fase yang berbeda, propagasi ketiga saraf

Dari hasil pada Gambar 52. dapat disimpulkan bahwa untuk tipe 2 berlaku kekuatan kopling εij merupakan faktor utama dalam fenomena sinkronisasi pada propagasi saraf tipe 1 maupun 2. 4.4.2 Solusi numerik pada 4 saraf terkopel Bahasan mengenai kopling 4 saraf ini adalah lanjutan dari sub bab 4.4.1 dengan tujuan agar lebih memahami fenomena kopling pada saraf. Pada sistem kopling n=4, maka persamaan (36) menjadi.

41 N . =− NX

. ∞ ( . )( .

−

)

− !. ( . − ) − "( . − ") . . + $z { sin (~X) + $ a . + $µ¹ ( / , 2 , 3 )(62. G) .( ) N!. !∞ . − !. = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (62. I) .( ) NX '( . N / / ) =− ∞ ( / )( / − NX − !/ ( / − ) − "( / − ") / / + $z { sin(~X) + $ a / ( + $µ¹ . , 2 , 3 )(62. P) /( ) N!/ !∞ / − !/ = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (62. N) /( ) NX '( / N 2 2 ) =− ∞ ( 2 )( 2 − NX − !2 ( 2 − ) − "( 2 − ") 2 2 + $z { sin(~X) + $ a 2 ( + $µ¹ . , / , 3 )(62. O) 2( ) N!2 !∞ 2 − !2 = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (62. @) 2( ) NX '( 2 N 3 3 ) =− ∞ ( 3 )( 3 − NX − !3 ( 3 − ) − "( 3 − ") 3 3 + $z { sin (~X) + $ a 3 + $µ¹ ( . , / , 2 )(62. ) 3( ) N!3 !∞ 3 − !3 = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ (62. ℎ) 3( ) NX '( 3

/ $µ¹ = −(b/

−

µ ) -ℎ/. q/. º

1

» 1 + expjσ(b. − t/ )m 1 + ℎ/2 q/2 º » 1 + expjσ(b2 − t/ )m 1 »´ (62. I) + ℎ/3 q/3 º 1 + expjσ(b3 − t/ )m 2 $µ¹ = −(b2 1 − µ ) -ℎ2. q2. º » 1 + expjσ(b. − t2 )m 1 + ℎ2/ q2/ º » 1 + expjσ(b/ − t2 )m 1 + ℎ23 q23 º »´ (62. P) 1 + expjσ(b/ − t2 )m 3 $µ¹ = −(b3 1 − µ ) -ℎ3. q3. º » 1 + expjσ(b. − t3 )m 1 + ℎ3/ q3/ º » 1 + expjσ(b/ − t3 )m 1 + ℎ32 q32 º »´ (62. N) 1 + expjσ(b2 − t3 )m Bentuk diagram pada kopling 4 saraf ini adalah. neuron 1

neuron 3

neuron 2

neuron 4

Dengan fungsi Isyn adalah.

. $µ¹ = −(b.

−

µ ) -ℎ./ q./ º

1

» 1 + expjσ(b/ − t. )m 1 + ℎ.2 q.2 º » 1 + expjσ(b2 − t. )m 1 + ℎ.3 q.3 º »´ (62. G) 1 + expjσ(b3 − t. )m

Gambar 53. Sistem 4 saraf terkopel. Hasil yang didapatkan pada simulasi propagasi tipe 1 dan 2 dengan asumsi bahwa tiap-tiap saraf pada sistem terhubung dengan kekuatan kopling yang seragam yaitu ε=0.5 ditampilkan pada Gambar 54.

42

neuron 1

4 Coupled Class 1 Excitability

neuron 4

neuron 3

neuron 2

coupled hij=1, epsij=0.5, not same phase 100

(a)

0 -100

0

100

200

full coupled hij=1, epsij=0.5, synchron

300

400

500

600

700

800

900

1000

50

(b)

0 -50

0

100

200

300

400

500

600

700

800

900

1000

300

400

500

600

700

800

900

1000

300

400

500 time (ms)

600

700

800

900

1000

4 Coupled Class 2 Excitability 100

coupled hij=1, epsij=0.5, not same phase

(c)

0 -100

0

100

200

dengan nilai ε13=ε31=0.5, untuk saraf 1 dan 2 serta 3 dan 4 terkopel namun memiliki kekuatan ikatan yang berbeda yaitu ε21=ε43=0.25 dan ε12=ε34=1.25, sedangkan untuk saraf 1 dan 4 serta 2 dan 3 tidak terkopel sempurna (ε41=ε23=0.5, ε14=ε32=0). Antara saraf 2 dan 4 tidak terkopel sama sekali (ε24=ε42=0).Hasil dari sistem pada Gambar 55. Adalah sebagai berikut.

full coupled hij=1, epsij=0.5, synchron 50 0 -50

(d) 0

100

200

Gambar 54. Sinkronisasi 4 saraf terkopel tipe 1 dan 2. (a) tipe 1 beda fase (b) tipe 1 sefase (c) tipe 2 beda fase (d) tipe 2 sefase. Keadaan sinkronisasi pada 4 saraf terkopel tidak jauh berbeda dengan kopel sebelumnya hanya yang membedakan adalah jumlah dari saraf yang terkopel. Pada tipe 1 maupun 2, sinkronisasi phase locking dapat terjadi dengan syarat nilai εij pada tiap-tiap ikatan adalah sama. Untuk dps terjadi saat keempat saraf memiliki fase awal propagasi yang berbeda yaitu dengan nilai awal potensial masing-masing untuk saraf 1 sampai 4 adalah -60 mV, -40 mV, -20 mV, dan 0 V. Untuk sps, fase awal propagasi keempat saraf disamakan yaitu memiliki nilai awal potensial 0 mV. Agar lebih memahami keadaan sinkronisasi pada 4 saraf terkopel ini dibangun suatu sistem diagram sebagai berikut. neuron 1

neuron 3

neuron 2

neuron 4

Gambar 55.Variasi sistem 4 saraf terkopel. Dari diagram tersebut dapat dilihat bahwa antara saraf 1 dan 3 terkopel

Gambar 56.Variasi propagasi 4 saraf terkopel (a) tipe 1 (b) tipe 2. Berdasarkan Gambar 56, hasil yang didapatkan adalah baik tipe 1 maupun 2 tidak mencapai sinkronisasi yang sempurna meskipun empat jenis saraf pada kedua tipe memiliki fase awal yang sama (Vinit = 0 mV). Pada tipe 1 dan 2, saraf 1 dan 3 hampir tersinkronisasi sempurna. Jika dilihat dari diagram fase nya, dengan garis merah (saraf 1) dan biru (saraf 3) saling berhimpitan. Ini juga terjadi pada saraf 2 (garis hijau) dan 4 (garis jingga) yang saling berhimpitan. Jika dilihat dari diagram sistem kopel pada Gambar 35., hasil sinkronisasi ini terjadi karena pada saraf 1 dan 3 yang saling berhimpitan memiliki ikatan yang kuat (hij=1) dengan nilai ε yang sama.Pada saraf 2 dan 4, meskipun tidak terkopel, tetapi memiliki fase awal yang sama dan nilai parameter saraf yang sama. Hal ini yang menyebabkan propagasi kedua saraf ini saling berhimpitan. Antara saraf 1 dan 2 serta saraf 3 dan 4 tidak saling berhimpitan walaupun terkopel. Hal ini terjadi karena

43

kedua saraf tersebut memiliki kekuatan kopel yang berbeda (εi≠εj). Telah dibahas keseluruhan simulasi model saraf Morris-Lecar dengan mengambil propagasi eksitasi neural properties tipe 1 dan 2. Jenis simulasi yang dilakukan pada tipe 1 dan 2 adalah dengan meninjau nilai parameter arus terapan Iapp DC tetap, DC bergantung waktu, dan AC. Untuk mengetahui makna kualitatif dari propagasi saraf ini dilakukan analisis sistem dinamik dengan mencari titik keseimbangan dan nilai eigennya. Pada akhir pembahasan, dilakukan simulasi sistem banyak saraf (n ≥2) dengan nilai arus terapan bergantung waktu AC.

V. SIMPULAN DAN SARAN 5.1 Simpulan Telah dilakukan analisis model saraf Morris-Lecar untuk propagasi tipe 1 dan 2 (class 1 and 2 excitability). Analisis pada propagasi tipe 1 dan 2 dilakukan dengan varaisi nilai arus terapan. Pada arus terapan DC konstan, tipe 1 dan 2 dapat memiliki propagasi yang stabil masing-masing saat nilai arus (Iapp=40 µA dan Iapp=50 µA) yaitu 50 µA dan 55 µA. Pada arus DC bergantung waktu, tipe 1 dan 2 memiliki batasan pita arus ekstasi yang berbeda. Tipe 2 memiliki daerah arus eksitasi yang lebih sempit (50-250 µA) dibandingkan tipe 1 (125-350 µA). Untuk pengaruh arus AC bergantung waktu, pada tipe 1 tidak mengalami fenomena burst melainkan hanya mengubah frekuensi eksitasi saja. Sedangkan pada tipe 2, dapat mengalami spike dan burst dengan variasi nilai ω. Analisis sistem dinamik pada tipe 1 dan 2 menghasilkan titik kritis saddle tidak stabil untuk tipe 1 dan focus stabil untuk tipe 2. Jenis bifurkasi pada Iapp tetap untuk tipe 1 adalah saddle-node on invariant circle (SNIC) sedangkan untuk tipe 2 adalah Andronov-Hopf. Pada nilai Iapp DC bergantung waktu, untuk tipe 1 mengalami perubahan bifurkasi dari unstable SNIC menjadi stable Andronov-Hopf. Sedangkan untuk tipe 2 hanya mengalami perubahan

kestabilan, dari unstable subcriticalAndronov-Hopf menjadi stable supercritical-Andronov-Hopf Saat nilai Iapp merupakan arus AC. Kedua tipe tidak mengalami perubahan jenis maupun kestabilan titik kritis,tetapi mengalami perubahan tingkat kestabilan. Pada tipe 1, nilai eigen menuju keadaan critical point (nilai eigen=0), sedangkan pada tipe 2 mengalami perubahan ke arah kestabilan yang semakin tinggi (nilai suku real eigen definit -∞). Pada arus AC ini pula diikiuti oleh pergeseran grafik garis nol untuk V (V nullcline), mengikuti osilasi dari arus sinusoidal yang menyebabkan letak titik keseimbangan ikut berosilasi. Model saraf terkopel yang dibangun merupakan suatu synaptic coupled dengan tipe propagasi inhibitory. Dari hasil simulasi untuk saraf n=2,3,4 jika dilihat dari skala waktu propagasinya didapatkan kesimpulan bahwa Jika dua saraf memiliki perbedaan fase propagasi, maka akan lebih mudah tersinkronisasi dibandingkan dengan dua saraf yang berbeda skala waktu propagasinya. Sedangkan untuk karakteristik kopling meliputi kekuatan ikatan dan keterhubungannya, untuk peristiwa sinkronisasi lebih dominan dipengaruhi oleh kekuatan kopling εij. Saraf dapat tersinkronisasi meskipun satu sama lain tidak saling terhubung (hi≠hji). Saat saraf terkopel (hi=hj=1) dan tersinkronisasi (εi=εj), keadaan ini disebut phase locking. Saat tercapai phase locking ada dua kemungkinan yang terjadi dilihat dari fase propagasi tiap saraf. Kopling saraf akan mengalami different phase synchronization apabila memiliki fase yang berbeda atau same phase synchronization apabila memiliki fase yang sama. 5.2 Saran Dalam penelitian ini masih terbatas pada beberapa analisis saja. Sebagai contoh untuk Iapp, arus DC bergantung waktu yang dipakai merupakan suatu fungsi linier, masih ada bentuk fungsi bergantung waktu lain seperti fungsi esponensial dan logaritmik. Fungsi AC bergantung waktu merupakan fungsi sinusoidal satu orde. Variasi

44

fungsi ini dapat dilakukan dengan fungsi sinusoidal 2 orde atau lebih atau bentuk fungsi AC lain seperti gelombang kotak, atau gergaji. Analisis sistem dinamik hanya terbatas pada variasi nilai parameter Iapp saja dan tidak meninjau variasi untuk parameter lain seperti gCa, gK, gL, dan konduktansi C membran. Selain itu, pada penelitian ini hanya menganalisis 2 bentuk propagasi saja yaitu tipe 1 dan 2. Sedangkan bentuk propagasi saraf memiliki banyak variasi seperti integrator, resonator, bistability, spike latency, periodic bursting, dan lainnya. Pada saraf terkopel, analisisnya hanya dibatasi pada karaktersitik sinkronisasi itu sendiri (kekuatan kopling, keterhubungan, fase propagasi, dan tidak membahas karaktersitik dinamik pada sistem terkopel utnuk banyak saraf. untuk kasus yang lebih kompleks, dapat dilakukan analisis propagasi saraf bergantung suhu baik single cell maupun coupled cell serta karakteristik dinamiknya.

5

6

7

8

9

10

11

DAFTAR PUSTAKA 1

2

3

4

Feriyawati L. Anatomi Sistem Saraf dan Peranannya dalam Regulasi Kontraksi Otot Rangka. [Disertasi]. Medan: Fakultas Kedokteran Universitas Sumatra Utara;2006. Keneer J, Sneyd J. Antman SS, Marsden JE, Sirovich L, editor. .Interdisciplinary Applied Mathematics, Mathematical Biology, Mathematical Physiology, System Physiology. Second ed. USA: ©Springer Science+Business Media, LLC.2008. Izhikevich EM. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Massachusetts:The MIT Press;2007. Corson N, Aziz-Zlaoui MA. Dynamics and complexity of Hindmarsh-Rose neuronal systems. France : Laboratoire de

12

13

14

15

Mathématques Appliquées du Havre. Le Havre.2008 Izhikevich EM.. Simple model of spiking neurons.IEEE Transactions on neural network,in press 2003;vol. 14:6 Hodgkin A L, Huxley A F. A Quantitative description of membrane current and application to conduction and excitation in nerve 1954:Journal Pysiol, 177 pp. Alatas H. Mathematics of Physics Lecturer. Bogor:Departemen Fisika Fakultas MIPA IPB;2008 Boas ML. Mathematical Methods in The Physical Science.New York:John wiley and Sons. Inc;1983 Izhikevich EM, Desai NS, Walcott EC, Hoppensteadt FC. Burst as a unit of neural information: Selective Communication via Resonance Trends in Neurosco., 2003;vol. 26, pp:161-167. Buzsaki G.. Rhythms of the brain, Oxford University Press.2006:vol. 32:21-25 Batista CAS et al. Global synchronization of bursting neurons in clustered networks. Santo Andr´e, S˜ao Paulo, Brazil.2002 Ueta T, Chen G). On synchronization and control of coupled Wilson Cowan neural oscillators. International Journal of Bifurcation and Chaos. World Scientic Publishing Company.2001;vol. 13, No. 1 163175. Izhikevich EM. Weakly pulsecoupled oscillators, fm interactions, synchronization , and oscillatory associative memory. IEEE transactions on neural networks .1999;vol. 10:3 Kavasseri RG, Ngarajan R. Synchronization in electrically coupled neural Networks.USA.2003 Agnon-Snir H, Carr CE, Rinzel J. The role of dendrites in auditory

45

16.

17

18

19

20

21

22

23

24

25

coincidence detectors. Nature.1998;393: 268 - 272. Kitajima H, Kurths J. Forced synchronization in Morris–Lecar neuron. International Journal of Bifurcation and Chaos.1997;Vol. 17:10 3523–3528. Kavasseri RG.Bifurcations in a synaptically coupled Morris-Lecar neuron model. Department of Electrical and Computer Engineering North Dakota State University, Fargo,2003:ND 58105 – 5285. Hoppensteadt FC, Izhikevich EM. Weakly Connected neural networks. Springer-Verlag 2007. Izhikevich EM. Class 1 neural excitability, conventional synapses, weakly connected networks, and mathematical foundations of pulsecoupled models, IEEE Trans. Neural Networks.1999;10:499-507. David DC. Introduction to bifurcation theory. The university of Texas at Austin. Institute for fusion study;1997 Sato YD , Shiino M, Spiking neuron models with axcitatory or inhibitory couplings and synchronization phenomena, Physical Review 2002;vol. 66:041903, pp.1-13. Izhikevich EM, Kuramoto Y. Weakly coupled oscillators. Elsevier Ltd. All rights reserved.2006 Izhikevich EM, Hoppensteadt FC. Slowly coupled oscillators: phase dynamics and synchronization. Society for Industrial and Applied Mathematics.2003;Vol. 63: 6, pp. 1935–1953. Christopher PF et al. Interdisciplinary Applied Mathematics, Mathematical Biology: Computational Cell Biology. Joel Edward Kaiser. Copyright of Material;1996 Izhikevich EM.Which Model to use for cortical spiking neuron?. IEEE Trans on Neural Network, IN PRESS, SUMMER.2004

26 Kupfermann I. Functional studies of cotransmission. Physiol.1991:Rev. 71: 683. 27 Hoppensteadt FC, Izhikevich EM. Oscillatory neurocomputers with dynamic connectivity.Physical review letters.1999;82:14. 28 Arecchi FT. Chaotic neuron dynamics, Synchronization,and Feature Binding: Quantum Aspects.Mind and Matter.2005;1(1):15-43. 29 Izhikevich EM. Multiple cusp bifurcation.Neural Network.1998;11:495-508. 30 Oh M, Matveev V. Loss of phase locking in non-weakly coupled inhibitory networks of type-I model neurons. comput Neurosci. 2008: DOI 10.1007/s10827-0080112-8. 31 Izhikevich EM. Class 1 neural excitability, conventional synapses, weakly connected network, and mathematical foundations of pulsecoupled models. IEEE Trans on Neural Network.1999;10:3. 32 Belykh I, de Lange E, Hasler M. Synchronization of bursting neuron: what matters in the network topology.physical review letters.2005;94:188101 33 Arrecchi TF. Chaotic neuron dynamics, synchronization, and feature binding: quantum aspect. Mind and Matter;2006:vol 1(1), pp. 15-43 34 Rabunal JR, Dorado J. Artificial Neural Networks in Real-Life Appication. IDEA group publishing.2006

LAMPIRAN

47

Lampiran 1. Diagram alir penelitian

Mulai

Pencaraian dan Penelusuran Pustaka

Pustaka Lengkap dan Siap Penelitian

Analisis Numerik Sistem Dinamik PD Izhikevich

Variasi Nilai Parameter Iapptetap, DC dan AC BergantungWaktu

Analisis Numerik Sistem PD Morris-Lecar dengan RK-4

Pembahasan dan Analisis Lanjut Analisis Numerik untuk Nilai danVektor Eigen Analisis Numerik Sinkronisasi Model Saraf ML Terkopel (n=1,2,3)

Penyususnan Laporan Penelitian

Ya Tujuan Penelitian Tercapai

Tidak

SELESAI

48

Lampiran 2. Metode Rungge-Kutta orde 4 (RK-4) Jika terdapat suatu persamaan differensial biasa (PDB) dengan variabel perubahan x terhadap waktu t. =

,

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 1

Dengan mendiskritkan persamaan 1 saat limit dt→0 maka. =

lim

→0

,

+1 −

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 2 =

,

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 3

Dengan mengganggap bahwa dt merupakan suatu increament h.maka, +1 =

+ℎ

,

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 4

Untuk metode RK-4, dilakukan ketelitian hingga orde 4 pada deret Taylor. Dengan menkoreksi penambahan fungsihf(x,t).yaitu sebagai berikut. 1

=ℎ

2

=ℎ

3

=ℎ

4

=ℎ

,

1 + ℎ, 4 3 + ℎ, 8 +

12 ℎ, 13

+ +

1 4

3 32

+

1 1

+

1932 2197

9 32 1

−

2

7200 2197

2

+

7296 2197

3

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 5

Sehingga persamaan (4) menjadi. +1 =

+

1+2 2+2 3+ 4 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 6 6

Dengan nilai s adalh increament yang lebih teliti dari h: *

*

∈ ℎ + ∈ ℎ + !=" ) = 0,840896 " ) ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 7 | &'( | 2| &'( | dimana∈ adalah toleransi nilai error.

49

Lampiran 3 Metode Rungge-Kutta orde 4 untuk model Morris-Lecar Model Morris-Lecar (1981) dengan dua variabel dimensional V dan W adalah sebagai berikut. ,

-

4

= −./0 12 - - − -/0 − .3 4 - − -3 − .5 - − -5 + 6077 ∙∙∙∙∙∙∙∙∙∙ 1

=

4∞ - − 4 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 2 9: -

Dengan menempatkan laju perubahan dV dan dW masing-masing di ruas kiri, maka. -=

;, :

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 3

4 = . ;, :

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 4

Dengan nilai f(v,w) dan g(v,w )masing-masing:

−./0 12 - - − -/0 − .3 4 - − -3 − .5 - − -5 + 6077 ∙∙∙∙∙∙∙ 5 , 42 - − 4 . ;, : = ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 6 9< ;, : =

Persamaan (3) dan (4) akan didiskritisasi sehingga menjadi. -

4

+1 = -

+ ℎ -, 4,

+1 =

+ ℎ -, 4,

∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 7 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 8

Model ML dengan metode RK-4 menjadi. Untuk V: 1

=ℎ

2

=ℎ =

3

=ℎ

4

=ℎ

-

,-

1

+

+1 = -

2

=ℎ

=ℎ =

1>

3 3 + ℎ, - + 8 32

,4

1

+

12 1932 ℎ, - + 13 2197 +

Dan untuk W: 1

1

+ 4 ℎ, - + 4

1

9 32 1

−

2

7200 2197

2

+

7296 2197

3

1+2 2+2 3+ 4 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 9 6

1

+ 4 ℎ, 4 + 4

1>

50

3

=ℎ

4

=ℎ

4

3 3 + ℎ, 4 + 8 32 +

1

+

12 1932 ℎ, 4 + 13 2197

+1 =4

+

9 32 1

−

2

7200 2197

2

+

7296 2197

3

1+2 2+2 3+ 4 ∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙ 10 6

Kemudian metode RK-4 untuk model ML ini akan disimulasikan dengan menggunakan MATLAB untuk mencari solusi secara numerik.

51

Lampiran 4. Program untuk mencari solusi numerik propagasi saraf function [t,V,w]=class1currenttime(mrva,mrwa,dt,n) t =[]; V =[]; w =[]; V(1)=0; %initial value of membrane Potential w(1)=0; %initial value of recovery variable tau = dt; %timespan (cacahan) t(1)= 0; for k=1:n t(k+1)=t(k)+tau; k1v=tau*feval(mrva, t(k), V(k), w(k)); k1w=tau*feval(mrwa, t(k), V(k), w(k)); k2v=tau*feval(mrva, t(k)+tau/2, V(k)+k1v/2, w(k)+k1w/2); k2w=tau*feval(mrwa, t(k)+tau/2, V(k)+k1v/2, w(k)+k1w/2); k3v=tau*feval(mrva, t(k)+tau/2, V(k)+k2v/2, w(k)+k2w/2); k3w=tau*feval(mrwa, t(k)+tau/2, V(k)+k2v/2, w(k)+k2w/2); k4v=tau*feval(mrva, t(k)+tau, V(k)+k3v/2, w(k)+k3w); k4w=tau*feval(mrwa, t(k)+tau, V(k)+k3v/2, w(k)+k3w); V(k+1)=V(k)+(k1v+ 2*k2v+ 2*k3v+ k4v)/6; w(k+1)=w(k)+(k1w+ 2*k2w+ 2*k3w+ k4w)/6; end subplot(2,2,1); plot(t,V1,'-r'); hold on; plot(t,V2,'-g'); hold on; plot(t,V3,'-b'); hold on; plot(t,V4,'-y'); %tau=[0:1:1800]; %alpha=0.056; %Im=5; %Im=10; %Iapp=Im*alpha*tau; %Iapp=Im*sin(alpha*tau); %plot(t,V2,'-g'); %figure subplot(2,2,2) plot (V1,w1,'-r');

hold on; plot (V2,w2,'-g'); hold on plot (V3,w3,'-b'); hold on; plot (V4,w4,'-y'); subplot(2,2,4) plot3 (V1,w1,t,'-r'); hold on; plot3 (V2,w2,t,'-g'); hold on; plot3 (V3,w3,t,'-b'); hold on; plot3 (V4,w4,t,'-y'); subplot(2,2,3); tau=[0:1:1000]; alpha=0.016; %Im=8; Im=10; %Iapp=Im*alpha*tau; Iapp=Im*sin(alpha*tau)+55; plot(tau,Iapp); %CopyRight of Adam Sukma Putra %[t,V,w]=class1currenttime(@mrva,@mrwa,0. 01,1000) %Pemanggilanfungsi parameter V functionmrva=mrva(t,V,w) C=20; %membrane capacitance (mikroFarad/cm^2) %Iapp=-14; %applied dc contant current (microAmpere/cm^2) I=50; %class I I=40 mikroVolt %I=55; %class II I=50 mikroVolt

%I for once spike %I for once spike

%omega=0; %no AC current omega=0.051; %class I %omega=0.016; %class II Im=50; %class I %Im=10; %class II Iapp=Im*sin(omega*t)+I; %injection AC current %Iapp=Im*omega*t; %injection DC current %significant parameter Va=-1.2; %in mV Vb=18; %in mV

52

%reamaining parameter VCa=120; %initial value of Ca channel potential (mV) VK=-80; %initial value of K channel potential (mV) VL=-60; %initial value of Leakage channel potential (mV) gCa=4; %initial value of CaCahnnel conductance (mS/cm^2) gK=8; %initial value of K channel conductance (mS/cm^2) gL=2; %initial value of Leak channel conductance (mS/cm^2) %remaining variable ICa=gCa*(0.5*(1+tanh((V-Va)/Vb)))*(VVCa); IK=gK*w*(V-VK); IL=gL*(V-VL); mrva=(-ICa-IK-IL+Iapp)/C; %Pemanggilanfungsi parameter W

functionmrwa=mrwa(t,V,w) psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point Vc=12; %Class I %Vc=2; %Class II %Vc=18; %Resonator %Vc=-3; %Integrator mrwa=(0.5*(1+tanh((V-Vc)/Vd))w)*(psi)*cosh((V-Vc)/(2*Vd));

Lampiran 5. Program untuk mencari nilai akar keseimbangan %pencarianakarsistem Morris-Lecar %using solve syntax %class 1 %[V,w] = solve('(0.5*(1+tanh((V-12)/17.4))w)*(1/15)*cosh((V-12)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+50)/20=0')

%class 2 %[V,w] = solve('(0.5*(1+tanh((V-2)/17.4))w)*(1/15)*cosh((V-2)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+55)/20=0') %class 1 Vc=18 %[V,w] = solve('(0.5*(1+tanh((V-18)/17.4))w)*(1/15)*cosh((V-18)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+50)/20=0') %class 1 Vc=-3 %[V,w] = solve('(0.5*(1+tanh((V+3)/17.4))w)*(1/15)*cosh((V+3)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+50)/20=0') %class 2 Vc=18 %[V,w] = solve('(0.5*(1+tanh((V-18)/17.4))w)*(1/15)*cosh((V-18)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+55)/20=0') %class 2 Vc=-3 %[V,w] = solve('(0.5*(1+tanh((V+3)/17.4))w)*(1/15)*cosh((V+3)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+55)/20=0') %class 1 DC current time %[V,w] = solve('(0.5*(1+tanh((V-12)/17.4))w)*(1/15)*cosh((V-12)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)2*(V+60)+(5*0.011*2300))/20=0') %class 2 DC current time %[V,w] = solve('(0.5*(1+tanh((V-2)/17.4))w)*(1/15)*cosh((V-2)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)2*(V+60)+(10*0.016*1800))/20=0') %class 1 AC current time %[V,w] = solve('(0.5*(1+tanh((V-12)/17.4))w)*(1/15)*cosh((V-12)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)2*(V+60)+(8*sin(0.011*425)+50))/20=0') %class 2 AC current time %[V,w] = solve('(0.5*(1+tanh((V-2)/17.4))w)*(1/15)*cosh((V-%2)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)%2*(V+60)+(10*sin(0.016*300)+55))/20=0') %result the roots of model class 1

53

%V=-5.375751747601989319002873223103 %w=0.1194958661895999719291973003847

%Iapp=-14; %applied dc contant current (microAmpere/cm^2)

%result the roots of model class 2 %V=-5.7730477497112001045947732959849 %w=0.11477389337488650867727252540234

%I=50; %class I I=40 mikroVolt I=55; %class II I=50 mikroVolt

%result the roots of model class 1 vc=18 %V=-7.7688412951999231678106681234789 %w=0.04917541536377540186546723826383 7 %result the roots of model class 1 vc=-3 %V=-37.944399395986573083934449704091 %w=0.01769538423492359794552919543552 2 %result the roots of model class 2 vc=18 %V=-8.2366207736868930934453826856985 %w=0.04672146056708466811143485370889 5 %result the roots of model class 2 vc=-3 %V=-36.216171408244318706082768732896 %w=0.02150036006662151829129503110092 8

%I for once spike %I for once spike

%t=700; %t=1500; %t=2000; %t=2300;

%class 1 DC current time

%t=300; %t=1000; %t=1600; %t=1800;

%Class 2 DC current time

%t=150; %t=280; %t=425;

%Class 1 AC current time

t=100; t=200; t=300;

%Class 2 AC current time

%class 1 DC current time dependent t=700 %V0=-4.529591984400099139054741692399; %w0=0.1301137824219026230818754922359 3; %class 1 DC current time dependent t=1500 state (B) %V0=8.7000878923783390351159092599439; %w0=0.0847632694996214267651510264315 13;

%omega=0; %no AC/DC current time %omega=0.011; %class I omega=0.0016; %class II %Im=8; %class I Im=10; %class II Iapp=Im*sin(omega*t)+I; %injection AC current %Iapp=Im*omega*t; %injection DC current %significant parameter

Lampiran 6. Program untuk mencari matriks jacobian dan nilai eigen

Va=-1.2; %in mV Vb=18; %in mV

% 'pencarian matrix jacobiandannilaieigen' %morris-lecar %fungsi diff untukmencariturunanfungsi %fungsi solve untukmencariakar %fungsidsolveuntukmencarisolusi PDB

%reamaining parameter VCa=120; %initial value of Ca channel potential (mV) VK=-80; %initial value of K channel potential (mV) VL=-60; %initial value of Leakage channel potential (mV)

%parameter V C=20; %membrane capacitance (mikroFarad/cm^2)

gCa=4; %initial value of CaCahnnel conductance (mS/cm^2)

54

gK=8; %initial value of K channel conductance (mS/cm^2) gL=2; %initial value of Leak channel conductance (mS/cm^2) %parameter w psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point %Vc=12; %Class I %node stable twin negative eigen value Vc=2; %Class II %node unstable twin positive eigen value %Vc=18; %integrator %saddle (I) saddle (II) paired eigen value %Vc=-3; %resonator %imajiner (I) imajiner (II) eigen value syms V w d1=((-gCa*(0.5*(1+tanh((V-Va)/Vb)))*(VVCa))-(gK*w*(V-VK))-(gL*(V-VL))+Iapp)/C; d2=(0.5*(1+tanh((V-Vc)/Vd))w)*(psi)*cosh((V-Vc)/(2*Vd));

b11=diff(d1,V); b12=diff(d1,w); %MatriksJacobian b21=diff(d2,V); b22=diff(d2,w); %mencariakar-akarsaatkesetimbangan %class 1 %[V,w] = solve('(0.5*(1+tanh((V-12)/17.4))w)*(1/15)*cosh((V-12)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+50)/20=0') %class 2 %[V,w] = solve('(0.5*(1+tanh((V-2)/17.4))w)*(1/15)*cosh((V-2)/(2*17.4))=0','(4*(0.5*(1+tanh((V+1.2)/8)))*(V-120)8*w*(V+80)-2*(V+60)+55)/20=0')

%syms x y z clear %Class 1 %V0=-5.375751747601989319002873223103; %w0=0.1194958661895999719291973003847; %Class 2 %V0=34.968338919656846612563231095693; %w0=0.0140742472565195467774536905483 16;

%class 1 Vc=-3 %V0=37.944399395986573083934449704091; %w0=0.0176953842349235979455291954355 22; %class 1 Vc=18 %V0=7.7688412951999231678106681234789; %w0=0.0491754153637754018654672382638 37; %class 2 Vc=-3 %V0=36.216171408244318706082768732896; %w0=0.0215003600666215182912950311009 28; %class 2 Vc=18 %V0=8.2366207736868930934453826856985; %w0=0.0467214605670846681114348537088 95; %class 1 DC current time dependent t=700 state (A) %Eigen value %V0=-4.529591984400099139054741692399; %0.3691 %w0=0.1301137824219026230818754922359 3; %-0.0085 %class 1 DC current time dependent t=1500 state (B) %V0=8.7000878923783390351159092599439; %0.3642 %w0=0.0847632694996214267651510264315 13; %-0.0336 %class 1 DC current time dependent t=2000 trantition state (T) %V0=12.529163879224132504050789578285; %-0.0933 + 0.2648i %w0=0.5152011725049665591408476684716 4; %-0.0933 - 0.2648i %class 1 DC current time dependent t=2300 state (C) %V0=13.11172960546653632198835035129; %-0.1057 + 0.2638i %w0=0.5319028530085752428063176633422 1; %-0.1057 - 0.2638i %class 2 DC current time dependent t=300 state (A) %V0=37.709692307111599436580454320388; %-0.0820 + 0.0345i

55

%w0=0.0103094389711560169121676653449 34; %-0.0820 - 0.0345i %class 2 DC current time dependent t=1000 state (B) %V0=5.9741094159465852443484849915905; %0.1480 + 0.0258i %w0=0.2856584586180664142411893091332 9; %0.1480 + 0.0258i %class 2 DC current time dependent t=1600 Trantition state (T) %V0=10.704777534439395884400770671172; %-0.1079 + 0.2342i %w0=0.7311665327781792089989664452585 5; %-0.1079 - 0.2342i %class 2 DC current time dependent t=1800 state (C) %V0=12.072181990203718746824827951892; %-0.1366 + 0.2214i %w0=0.7609185563276897463969430853177; %-0.1366 - 0.2214i %class 1 AC current time dependent t=150 state max current %V0=6.0208217973534591376859005443927; %0.3770 %w0=0.1119119166858133921084642939847 8; %-0.0192 %class 1 AC current time dependent t=280 state med current %V0=5.4139361590118737377051019599215; %0.3753 %w0=0.1190348398379777661399330917893 6; %-0.0152 %class 1 AC current time dependent t=425 state min current %V0=4.7794356681838288180890854950365; %0.3713 %w0=0.1268977791737155389580968417874; %-0.0105 %class 2 AC current time dependent t=100 state max current %V0=31.424598212411997057733474783538; %-0.0528 + 0.0487i %w0=0.0210020840498290323961759956021 65; %-0.0528 - 0.0487i %class 2 AC current time dependent t=200 state med current

%V0=35.188976592640549830511494321184; %-0.0715 + 0.0407i %w0=0.0137266410379793381147123107882 37; %-0.0715 - 0.0407i %class 2 AC current time dependent t=300 state min current %V0=38.929452806492129651055164715646; %-0.0866 + 0.0313i %w0=0.0089728835581524423292995258889 34; %-0.0866 - 0.0313i J=[b11 b12;b21 b22] %matriksjacobiandalambentukpersamaan Jacob=subs(J,{V,w},{V0,w0}) %matriksjacobiandalambentuknilai Lambdas=eig(Jacob) %mencarinilaieigen

Lampiran 7. Program untuk grafik garis nol V nulclines dan W nulclines Vb=18; Vc=2; Vd=17.4; gCa=4; gK=8; gL=2; VCa=180; VK=-120; VL=-40; Vo = -60:0.1:40; Iapps=Im*sin(alpha*300)+55; Nullslow = (0.5)*(1+tanh((Vo-Vc)/Vd)); minf = 0.5*(1+tanh((Vo-Va)/Vb)); Nullfast = -(Iapps + gL*(VL-Vo)+ gCa*minf.*(VCa-Vo))./(gK*(VK-Vo)); plot(Vo, Nullfast,'c',Vo,Nullslow,'r'); %CopyRight of Adam Sukma Putra

Lampiran 8. Program fterkopel n=2,3,4.

untuk

sara

Untuk n=2 Function[t,V1,w1,V2,w2]=class12coupled(mrv a1,mrwa1,mrva2,mrwa2,dt,n)

56

t =[]; V1 =[]; w1 =[]; V2=[]; w2=[]; V1(1)=-40; V2(1)=0; %initial value of membrane Potential w1(1)=0; w2(1)=0; %initial value of recovery variable tau = dt; %timespan (cacahan) t(1)= 0; for k=1:n t(k+1)=t(k)+tau; k1v1=tau*feval(mrva1, t(k), V1(k), w1(k), V2(k)); k1v2=tau*feval(mrva2, t(k), V2(k), w2(k), V1(k)); k1w1=tau*feval(mrwa1, t(k), V1(k), w1(k)); k1w2=tau*feval(mrwa2, t(k), V2(k), w2(k)); k2v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2, V2(k)+k1v1/2); k2v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2, V1(k)+k1v2/2); k2w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2); k2w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2); k3v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2, V2(k)+k2v1/2); k3v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2, V1(k)+k2v2/2); k3w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2); k3w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2); k4v1=tau*feval(mrva1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1, V2(k)+k3v1/2); k4v2=tau*feval(mrva2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2, V1(k)+k3v2/2); k4w1=tau*feval(mrwa1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1); k4w2=tau*feval(mrwa2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2); V1(k+1)=V1(k)+(k1v1+ 2*k2v1+ 2*k3v1+ k4v1)/6; V2(k+1)=V2(k)+(k1v2+ 2*k2v2+ 2*k3v2+ k4v2)/6; w1(k+1)=w1(k)+(k1w1+ 2*k2w1+ 2*k3w1+ k4w1)/6; w2(k+1)=w2(k)+(k1w2+ 2*k2w2+ 2*k3w2+ k4w2)/6; end

Untuk n=3

function [t,V1,w1,V2,w2,V3,w3]=classcoupled3(mrva1, mrwa1,mrva2,mrwa2,mrva3,mrwa3,dt,n) t =[]; V1 =[]; w1 =[]; V2=[]; w2=[]; V3=[]; w3=[]; V1(1)=0; V2(1)=0; V3(1)=0; %initial value of membrane Potential w1(1)=0; w2(1)=0; w3(1)=0; %initial value of recovery variable tau = dt; %timespan (cacahan) t(1)= 0; for k=1:n t(k+1)=t(k)+tau; k1v1=tau*feval(mrva1, t(k), V1(k), w1(k), V2(k), V3(k)); k1v2=tau*feval(mrva2, t(k), V2(k), w2(k), V1(k), V3(k)); k1v3=tau*feval(mrva3, t(k), V3(k), w3(k), V1(k), V2(k)); k1w1=tau*feval(mrwa1, t(k), V1(k), w1(k)); k1w2=tau*feval(mrwa2, t(k), V2(k), w2(k)); k1w3=tau*feval(mrwa3, t(k), V3(k), w3(k)); k2v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2, V2(k)+k1v1/2, V3(k)+k1v1/2); k2v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2, V1(k)+k1v2/2, V3(k)+k1v2/2); k2v3=tau*feval(mrva2, t(k)+tau/2, V3(k)+k1v3/2, w3(k)+k1w3/2, V1(k)+k1v3/2, V2(k)+k1v2/2); k2w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2); k2w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2); k2w3=tau*feval(mrwa3, t(k)+tau/2, V3(k)+k1v3/2, w3(k)+k1w3/2); k3v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2, V2(k)+k2v1/2, V3(k)+k2v1/2); k3v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2, V1(k)+k2v2/2, V3(k)+k2v2/2); k3v3=tau*feval(mrva3, t(k)+tau/2, V3(k)+k2v3/2, w3(k)+k2w3/2, V1(k)+k2v3/2, V2(k)+k2v3/2); k3w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2); k3w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2);

57

k3w3=tau*feval(mrwa3, t(k)+tau/2, V3(k)+k2v3/2, w3(k)+k2w3/2); k4v1=tau*feval(mrva1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1, V2(k)+k3v1/2, V3(k)+k3v1/2); k4v2=tau*feval(mrva2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2, V1(k)+k3v2/2, V3(k)+k3v2/2); k4v3=tau*feval(mrva3, t(k)+tau, V3(k)+k3v3/2, w3(k)+k3w3, V1(k)+k3v3/2, V2(k)+k3v3/2); k4w1=tau*feval(mrwa1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1); k4w2=tau*feval(mrwa2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2); k4w3=tau*feval(mrwa3, t(k)+tau, V3(k)+k3v3/2, w3(k)+k3w3); V1(k+1)=V1(k)+(k1v1+ 2*k2v1+ 2*k3v1+ k4v1)/6; V2(k+1)=V2(k)+(k1v2+ 2*k2v2+ 2*k3v2+ k4v2)/6; V3(k+1)=V3(k)+(k1v3+ 2*k2v3+ 2*k3v3+ k4v3)/6; w1(k+1)=w1(k)+(k1w1+ 2*k2w1+ 2*k3w1+ k4w1)/6; w2(k+1)=w2(k)+(k1w2+ 2*k2w2+ 2*k3w2+ k4w2)/6; w3(k+1)=w3(k)+(k1w3+ 2*k2w3+ 2*k3w3+ k4w3)/6; end

Untuk n=4 function [t,V1,w1,V2,w2,V3,w3,V4,w4]=classcoupled4( mrva1,mrwa1,mrva2,mrwa2,mrva3,mrwa3,mrv a4,mrwa4,dt,n) t =[]; V1 =[]; w1 =[]; V2=[]; w2=[]; V3=[]; w3=[]; V4=[]; w4=[]; V1(1)=0; V2(1)=0; V3(1)=0; V4(1)=0;%initial value of membrane Potential w1(1)=0; w2(1)=0; w3(1)=0; w4(1)=0;%initial value of recovery variable tau = dt; %timespan (cacahan) t(1)= 0; for k=1:n t(k+1)=t(k)+tau; k1v1=tau*feval(mrva1, t(k), V1(k), w1(k), V2(k), V3(k), V4(k)); k1v2=tau*feval(mrva2, t(k), V2(k), w2(k), V1(k), V3(k), V4(k));

k1v3=tau*feval(mrva3, t(k), V3(k), w3(k), V1(k), V2(k), V4(k)); k1v4=tau*feval(mrva3, t(k), V4(k), w4(k), V1(k), V2(k), V3(k)); k1w1=tau*feval(mrwa1, t(k), V1(k), w1(k)); k1w2=tau*feval(mrwa2, t(k), V2(k), w2(k)); k1w3=tau*feval(mrwa3, t(k), V3(k), w3(k)); k1w4=tau*feval(mrwa4, t(k), V4(k), w4(k)); k2v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2, V2(k)+k1v1/2, V3(k)+k1v1/2, V4(k)+k1v1/2); k2v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2, V1(k)+k1v2/2, V3(k)+k1v2/2, V4(k)+k1v2/2); k2v3=tau*feval(mrva3, t(k)+tau/2, V3(k)+k1v3/2, w3(k)+k1w3/2, V1(k)+k1v3/2, V2(k)+k1v3/2, V4(k)+k1v3/2); k2v4=tau*feval(mrva4, t(k)+tau/2, V4(k)+k1v4/2, w4(k)+k1w4/2, V1(k)+k1v4/2, V2(k)+k1v4/2, V3(k)+k1v4/2); k2w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k1v1/2, w1(k)+k1w1/2); k2w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k1v2/2, w2(k)+k1w2/2); k2w3=tau*feval(mrwa3, t(k)+tau/2, V3(k)+k1v3/2, w3(k)+k1w3/2); k2w4=tau*feval(mrwa4, t(k)+tau/2, V4(k)+k1v4/2, w4(k)+k1w4/2); k3v1=tau*feval(mrva1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2, V2(k)+k2v1/2, V3(k)+k2v1/2, V4(k)+k2v1/2); k3v2=tau*feval(mrva2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2, V1(k)+k2v2/2, V3(k)+k2v2/2, V4(k)+k2v2/2); k3v3=tau*feval(mrva3, t(k)+tau/2, V3(k)+k2v3/2, w3(k)+k2w3/2, V1(k)+k2v3/2, V2(k)+k2v3/2, V4(k)+k2v3/2); k3v4=tau*feval(mrva4, t(k)+tau/2, V4(k)+k2v4/2, w4(k)+k2w4/2, V1(k)+k2v3/2, V2(k)+k2v3/2, V3(k)+k2v3/2); k3w1=tau*feval(mrwa1, t(k)+tau/2, V1(k)+k2v1/2, w1(k)+k2w1/2); k3w2=tau*feval(mrwa2, t(k)+tau/2, V2(k)+k2v2/2, w2(k)+k2w2/2); k3w3=tau*feval(mrwa3, t(k)+tau/2, V3(k)+k2v3/2, w3(k)+k2w3/2); k3w4=tau*feval(mrwa4, t(k)+tau/2, V4(k)+k2v4/2, w4(k)+k2w4/2); k4v1=tau*feval(mrva1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1, V2(k)+k3v1/2, V3(k)+k3v1/2, V4(k)+k3v1/2); k4v2=tau*feval(mrva2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2, V1(k)+k3v2/2, V3(k)+k3v2/2, V4(k)+k3v2/2);

58

k4v3=tau*feval(mrva3, t(k)+tau, V3(k)+k3v3/2, w3(k)+k3w3, V1(k)+k3v3/2, V2(k)+k3v3/2, V4(k)+k3v3/2); k4v4=tau*feval(mrva4, t(k)+tau, V4(k)+k3v4/2, w4(k)+k3w4, V1(k)+k3v3/2, V2(k)+k3v3/2, V3(k)+k3v3/2); k4w1=tau*feval(mrwa1, t(k)+tau, V1(k)+k3v1/2, w1(k)+k3w1); k4w2=tau*feval(mrwa2, t(k)+tau, V2(k)+k3v2/2, w2(k)+k3w2); k4w3=tau*feval(mrwa3, t(k)+tau, V3(k)+k3v3/2, w3(k)+k3w3); k4w4=tau*feval(mrwa4, t(k)+tau, V4(k)+k3v4/2, w4(k)+k3w4); V1(k+1)=V1(k)+(k1v1+ 2*k2v1+ 2*k3v1+ k4v1)/6; V2(k+1)=V2(k)+(k1v2+ 2*k2v2+ 2*k3v2+ k4v2)/6; V3(k+1)=V3(k)+(k1v3+ 2*k2v3+ 2*k3v3+ k4v3)/6; V4(k+1)=V4(k)+(k1v4+ 2*k2v4+ 2*k3v4+ k4v4)/6; w1(k+1)=w1(k)+(k1w1+ 2*k2w1+ 2*k3w1+ k4w1)/6; w2(k+1)=w2(k)+(k1w2+ 2*k2w2+ 2*k3w2+ k4w2)/6; w3(k+1)=w3(k)+(k1w3+ 2*k2w3+ 2*k3w3+ k4w3)/6; w4(k+1)=w4(k)+(k1w4+ 2*k2w4+ 2*k3w4+ k4w4)/6; end subplot(2,2,1); plot(t,V1,'-r'); hold on; plot(t,V2,'-g'); hold on; plot(t,V3,'-b'); hold on; plot(t,V4,'-y'); %tau=[0:1:1800]; %alpha=0.056; %Im=5; %Im=10; %Iapp=Im*alpha*tau; %Iapp=Im*sin(alpha*tau); %plot(t,V2,'-g'); %figure subplot(2,2,2) plot (V1,w1,'-r'); hold on; plot (V2,w2,'-g');

hold on plot (V3,w3,'-b'); hold on; plot (V4,w4,'-y'); subplot(2,2,4) plot3 (V1,w1,t,'-r'); hold on; plot3 (V2,w2,t,'-g'); hold on; plot3 (V3,w3,t,'-b'); hold on; plot3 (V4,w4,t,'-y'); subplot(2,2,3); tau=[0:1:1000]; alpha=0.016; %Im=8; Im=10; %Iapp=Im*alpha*tau; Iapp=Im*sin(alpha*tau)+55; plot(tau,Iapp); %CopyRight of Adam Sukma Putra %[t,V1,w1,V2,w2,V3,w3,V4,w4]=classcoupled 4(@mrva1,@mrwa1,@mrva2,@mrwa2,@mrva 3,@mrwa3,@mrva4,@mrwa4,dt,n)

Lampiran 9. Program untuk saraf 1,2,3, dan 4 Saraf 1 function mrva1=mrva1(t,V1,w1,V2,V3,V4) C=20; %membrane capacitance (mikroFarad/cm^2) %Iapp=-14; %applied dc contant current (microAmpere/cm^2) %I=50; %class I I=40 mikroVolt I=55; %class II I=50 mikroVolt

%I for once spike %I for once spike

%omega=0; %no AC current %omega=0.011; %class I omega=0.016; %class II %Im=8; %class I Im=10; %class II Iapp=Im*sin(omega*t)+I; %injection AC current %Iapp=Im*omega*t; %injection DC current %Synchronization parameter Vc=2;

59

g12=1.25; g13=0.5; g14=0.5; teta=-40; lamb=0.01; c12=1; c13=1; c14=0; %if V2>teta % c=1; %else % c=1; %end H12=1/(1+exp(lamb*(V2-teta))); H13=1/(1+exp(lamb*(V3-teta))); H14=1/(1+exp(lamb*(V4-teta))); Isyn1=-(V1Vc)*(g12*c12*H12+g13*c13*H13+g14*c14* H14); Va=-1.2; %in mV Vb=18; %in mV %reamaining parameter VCa=120; %initial value of Ca channel potential (mV) VK=-80; %initial value of K channel potential (mV) VL=-60; %initial value of Leakage channel potential (mV) gCa=4; %initial value of CaCahnnel conductance (mS/cm^2) gK=8; %initial value of K channel conductance (mS/cm^2) gL=2; %initial value of Leak channel conductance (mS/cm^2) %remaining variable ICa=gCa*(0.5*(1+tanh((V1-Va)/Vb)))*(V1VCa); IK=gK*w1*(V1-VK); IL=gL*(V1-VL); mrva1=(-ICa-IK-IL+Iapp+Isyn1)/C; function mrwa1=mrwa1(t,V1,w1) psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point

%Vc=12; %Class I Vc=2; %Class II %Vc=18; %Resonator %Vc=-3; %Integrator mrwa1=(0.5*(1+tanh((V1-Vc)/Vd))w1)*(psi)*cosh((V1-Vc)/(2*Vd));

Saraf 2 function mrva2=mrva2(t,V2,w2,V1,V3,V4) %Synchronization parameter Vc=2; g21=0.25; g23=0.5; g24=0.5; teta=-40; lamb=0.01; c21=1; c23=1; c24=0; %if V1>teta % c=1; %else % c=1; %end H21=1/(1+exp(lamb*(V1-teta))); H23=1/(1+exp(lamb*(V3-teta))); H24=1/(1+exp(lamb*(V4-teta))); Isyn2=-(V2Vc)*(g21*c21*H21+g23*c23*H23+g24*c24* H24); %remaining variable ICa=gCa*(0.5*(1+tanh((V2-Va)/Vb)))*(V2VCa); IK=gK*w2*(V2-VK); IL=gL*(V2-VL); mrva2=(-ICa-IK-IL+Iapp+Isyn2)/C; function mrwa2=mrwa2(t,V2,w2) psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point %Vc=12; %Class I Vc=2; %Class II %Vc=18; %Resonator %Vc=-3; %Integrator

60

mrwa2=(0.5*(1+tanh((V2-Vc)/Vd))w2)*(psi)*cosh((V2-Vc)/(2*Vd));

mrwa3=(0.5*(1+tanh((V3-Vc)/Vd))w3)*(psi)*cosh((V3-Vc)/(2*Vd));

Saraf 3

Saraf 4

function mrva3=mrva3(t,V3,w3,V1,V2,V4) %Synchronization parameter Vc=2; g31=0.5; g32=0.5; g34=1.25; teta=-40; lamb=0.01; c31=1; c32=0; c34=1; %if V2>teta % c=1; %else % c=1; %end

function mrva4=mrva4(t,V4,w4,V2,V3,V1) %Synchronization parameter Vc=2; g41=0.5; g42=0.5; g43=0.25; teta=-40; lamb=0.01; c41=1; c42=0; c43=1; %if V2>teta % c=1; %else % c=1; %end

H31=1/(1+exp(lamb*(V1-teta))); H32=1/(1+exp(lamb*(V2-teta))); H34=1/(1+exp(lamb*(V4-teta)));

H41=1/(1+exp(lamb*(V1-teta))); H42=1/(1+exp(lamb*(V2-teta))); H43=1/(1+exp(lamb*(V3-teta)));

Isyn3=-(V3Vc)*(g31*c31*H31+g32*c32*H32+g34*c34* H34);

Isyn4=-(V4Vc)*(g41*c41*H41+g42*c42*H42+g43*c43* H43); %remaining variable ICa=gCa*(0.5*(1+tanh((V4-Va)/Vb)))*(V4VCa); IK=gK*w4*(V4-VK); IL=gL*(V4-VL);

%remaining variable ICa=gCa*(0.5*(1+tanh((V3-Va)/Vb)))*(V3VCa); IK=gK*w3*(V3-VK); IL=gL*(V3-VL); mrva3=(-ICa-IK-IL+Iapp+Isyn3)/C; function mrwa3=mrwa3(t,V3,w3) psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point %Vc=12; %Class I Vc=2; %Class II %Vc=18; %Resonator %Vc=-3; %Integrator

mrva4=(-ICa-IK-IL+Iapp+Isyn4)/C; function mrwa4=mrwa4(t,V4,w4) psi=1/15; %in s^-1 Vd=17.4; %in mV %Vc=4.6; %critical point %Vc=12; %Class I Vc=2; %Class II %Vc=18; %Resonator %Vc=-3; %Integrator mrwa4=(0.5*(1+tanh((V4-Vc)/Vd))w4)*(psi)*cosh((V4-Vc)/(2*Vd));

NOMENKLATUR No. 1 2

Simbol V W

Satuan mV mV

3

M∞(V) dan W∞ (V)

-

4 5

τω(V) V1, dan V3

mV

6

V2

mV

7 8 9 10 11 12 13 14 15

V4 mV ø s-1 ICa,IK, dan IL µA gCa, gK, dan gL mS/cm2 Iapp µA C µF/cm2 Iinit µA Imax µA α s-1

16 17 18 19 20 21

ω ε gij σ θj Vs

s-1 mS/cm2 mV mV

22

hij

-

23 24 25 26

Isyn J λ L

µA

Keterangan Potensial membran sel saraf Parameter pemulihan (depolarisasi saraf) yang berkaitan dengan normalisasi ion K+ Parameter fungsi konduksi, fungsi probabilitas keadaan terbuka dan tertutupnya saluran ionik yang sebanding dengan distribusi boltzman Fungsi waktu terbuka saluranK+. Nilai tengah saat arus ionic Ca2+ dan K+ada pada keadaan setengah teraktivasi (half activated) Konstanta potensial yang bertanggung jawab kepada loncatan potensial saat aktivasi Faktor kemiringan dari laju aktivasi ion K+ Skala waktu proses pemulihan Arus ion Ca+ , K+, dan leak current Nilai konduktansi membran Arus listrik terapan pada membran Kapasitansi membran sel saraf Nilai arus saat t=0 Nilai penambahan (gradien) arus maksimum tiap detik Penambahan yang bertanggung jawab atas besar kecilnya laju perubahan arus DC bergantung waktu Merupakan nilaifrekuensi masukan pada sinyal arus AC Kekuatan kopling Fungsi kopling Pengaturan laju kekuatan koling Potensial ambang tiap sel saraf ke-j Potensial pembalik dengan jenis hubungan sinaptik inhibitory Fungsi Heaviside yang menentukan apakah antara kedua saraf terhubung atau tidak Fungsi arus sinaptik Jacobian matriks Nilai eigen Vektor eigen