SIMULASI PROPAGASI IMPULS PADA SEL SARAF TERKOPEL MENGGUNAKAN MODEL FITZHUGH-NAGUMO
SUBIYANTO
DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011
ABSTRAK Subiyanto. Simulasi Propagasi Impuls pada Sel Saraf Terkopel Menggunakan Model FitzHugh-Nagumo. Dibimbing oleh Dr. Agus Kartono dan Dr.Ir. Irzaman, M.Si. Skripsi ini membuat sebuah model matematika untuk sel saraf terkopel. Model ini merupakan modifikasi dari model FitzHugh-Nagumo. Hasil simulasi model ini dapat menjelaskan perilaku sel saraf dalam menghantarkan impuls. Hasil simulasi dua sel saraf terkopel menunjukan bahwa arus minimal agar impuls dapat direspon oleh tubuh yaitu 0,41 μA, untuk model tiga sel saraf 0,47 μA, model empat sel saraf 0,54 μA, model delapan sel saraf 0,81 μA, dan model delapan belas sel saraf 1,4 μA. Dari hasil simulasi sel saraf ini didapatkan bahwa setiap penambahan satu sel saraf akan menambah arus minimal yang dapat direspon oleh tubuh itu sekitar 0,06 μA. Oleh karena itu, dengan menganggap bahwa jumlah sel saraf dalam tubuh yang berperan dalam penghantaran impuls itu sekitar satu juta sel saraf maka dapat diprediksi arus minimal yang dapat direspon oleh tubuh dalam penjalaran impuls sekitar 0,06 A, sedangkan arus yang kurang dari itu tidak akan mendapat respon dari tubuh. Kata kunci : model FitzHugh-Nagumo, sel saraf terkopel, impuls, respon
PERNYATAAN Dengan ini saya menyatakan bahwa skripsi yang berjudul Simulasi Propagasi Impuls pada Sel Saraf Terkopel Menggunakan Model FitzHugh-Nagumo adalah benar-benar hasil karya saya sendiri di bawah bimbingan Dr. Agus Kartono dan Dr.Ir. Irzaman, M.Si, dan belum pernah digunakan sebagai karya ilmiah pada perguruan tinggi atau lembaga manapun. Sumber informasi yang berasal atau dikutip dari karya yang diterbitkan maupun tidak diterbitkan dari penulis lain telah disebutkan dalam teks dan dicantumkan dalam Daftar Pustaka di bagian akhir skripsi ini.
Bogor,
Maret 2011
Subiyanto
Judul
: Simulasi Propagasi Impuls pada Sel Saraf Terkopel Menggunakan Model FitzHugh-Nagumo
Nama : Subiyanto NRP
: G74070017
Menyetujui,
(Dr. Agus Kartono)
(Dr.Ir. Irzaman, M.Si)
Pembimbing I
Pembimbing II
Mengetahui,
(Dr.Ir. Irzaman, M.Si) Ketua Departemen Fisika
Tanggal Lulus :
SIMULASI PROPAGASI IMPULS PADA SEL SARAF TERKOPEL MENGGUNAKAN MODEL FITZHUGH-NAGUMO
SUBIYANTO G74070017
Skripsi Sebagai salah satu syarat untuk memperoleh gelar Sarjana Sains pada Fakultas Matematika dan Ilmu Pengetahuan Alam Institut Pertanian Bogor
DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM INSTITUT PERTANIAN BOGOR BOGOR 2011
RIWAYAT HIDUP Penulis dilahirkan di Kediri, 5 Mei 1989 dari pasangan Supadi dan Erlina. Penulis merupakan anak ke empat dari empat bersaudara. Penulis menyelesaikan pendidikan Taman Kanak-Kanak hingga Sekolah Menengah Atas di Kabupaten Kediri yaitu TK Dharma Wanita Siman II, SD Negeri Siman IV, SMP Negeri 1 Kepung, dan SMA Negeri 1Pare. Penulis
melanjutkan
pendidikannya
ke
jenjang
perkuliahan di Departemen Fisika Institut Pertanian Bogor melalui jalur USMI. Selama kuliah, penulis aktif menjadi asisten praktikum Fisika TPB tahun 2008-2011, asisten praktikum mata kuliah Gelombang tahun 2009 dan 2011, asisten praktikum mata kuliah Fisika Komputasi tahun 2011. Penulis juga aktif sebagai kepala divisi Instrumentasi dan Teknologi di Himpunan Mahasiswa Fisika IPB tahun 2008-2009.
KATA PENGANTAR Puji syukur kepada Allah Subhanahu wa Ta‟ala, atas segala rahmat, nikmat kesehatan, kekuatan dan karunia-Nya sehingga penelitian untuk tugas akhir yang berjudul “Simulasi Propagasi Impuls pada Sel Saraf Terkopel Menggunakan Model FitzHughNagumo” dapat diselesaikan. Tugas akhir ini disusun sebagai salah satu syarat kelulusan program sarjana di Departemen Fisika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Pertanian Bogor. Pada kesempatan ini, penulis juga ingin mengucapkan terima kasih kepada orang tua dan saudara penulis yang selalu memberikan nasehat, bantuan dan semangat kepada penulis. Kepada Bapak Agus Kartono dan Bapak Irzaman sebagai pembimbing skripsi yang selalu memberikan motivasi dan semangat untuk menyelesaikan penelitian ini dan serta memberikan sebagian waktu selama berada di kampus hanya untuk berdiskusi mengenai penyusunan tugas akhir ini. Kepada Bapak Akhiruddin Maddu sebagai pembimbing akademik dan penguji yang telah banyak memberikan masukan kepada penulis ucapkan terima kasih. Kepada Ibu Mersi sebagai dosen penguji yang telah banyak memberikan masukan. Kepada rekan-rekan fisika beserta civitas akademik fisika lainnya yang telah banyak membantu penulis selama ini, penulis juga ucapkan terima kasih. Penulis menyadari bahwa penelitian tugas akhir ini jauh dari sempurna, oleh karena itu, saran dan kritik yang bersifat membangun sangat diperlukan bagi penulis. Semoga tugas akhir ini bermanfaat bagi semuanya.
Bogor, 20 Maret 2011 Penulis
DAFTAR ISI halaman DAFTAR GAMBAR .........................................................................................
iv
DAFTAR LAMPIRAN ......................................................................................
v
BAB 1 PENDAHULUAN ................................................................................
1
1.1 Latar Belakang ...................................................................................
1
1.2 Tujuan Penelitian ...............................................................................
1
1.3 Perumusan Masalah ...........................................................................
1
1.4 Hipotesis ............................................................................................
1
BAB 2 TINJAUAN PUSTAKA........................................................................
1
2.1 Sistem Saraf .......................................................................................
1
2.2 Model Hodgkin-Huxley .....................................................................
2
2.3 Model Fitzhugh-Nagumo ...................................................................
3
2.4 Model Matematika Baru Sel Saraf Terkopel ......................................
3
2.5 Mekanisme Penghantar Impuls ..........................................................
4
2.5.1 Penghantaran impuls melalui sel saraf ..........................................
4
2.5.1 Penghantaran impuls melalui sinapsis ..........................................
4
2.6 Motode Ode45 ...................................................................................
4
BAB 3 METODE PENELITIAN ......................................................................
5
3.1 Waktu dan Tempat Penelitian ............................................................
5
3.2 Peralatan ............................................................................................
6
3.3 Metode Penelitian ..............................................................................
6
3.3.1 Studi pustaka ................................................................................
6
3.3.2 Analisa numerik ...........................................................................
6
3.3.3 Analisa sifat fisis model yang telah dibuat ...................................
6
BAB 4 HASIL DAN PEMBAHASAN .............................................................
6
4.1 Simulasi Satu Sel Saraf ......................................................................
6
4.2 Model Matematika Sel Saraf Terkopel...............................................
7
4.2.1 Model dua sel saraf terkopel .........................................................
7
4.2.2 Model tiga sel saraf terkopel ........................................................
8
4.2.3 Model empat sel saraf terkopel .....................................................
8
4.2.4 Model delapan sel saraf terkopel ..................................................
9
4.2.4 Model enam belas sel saraf terkopel .............................................
10
4.3 Hasil Simulasi Model Sel Saraf Terkopel ..........................................
12
4.3.1 Hasil Simulasi dua Sel Saraf Terkopel .........................................
12
iii
4.3.2 Hasil Simulasi tiga Sel Saraf Terkopel .........................................
12
4.3.3 Hasil Simulasi empat Sel Saraf Terkopel .....................................
13
4.3.4 Hasil Simulasi delapan Sel Saraf Terkopel ...................................
13
4.3.5 Hasil Simulasi enam belas Sel Saraf Terkopel .............................
14
4.4 Analisa Model Sel Saraf Terkopel .....................................................
14
BAB 5 KESIMPULAN DAN SARAN .............................................................
15
5.1 Kesimpulan ........................................................................................
15
5.2 Saran ..................................................................................................
15
DAFTAR PUSTAKA ........................................................................................
15
LAMPIRAN.......................................................................................................
17
DAFTAR GAMBAR halaman Gambar 1. Bagian sel saraf ........................................................................................... 2 Gambar 2. Skema potensial aksi ................................................................................... 2 Gambar 3. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0 μA ..... 6 Gambar 4. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,2μA ... 6 Gambar 5. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,3μA ... 7 Gambar 6. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,35μA . 7 Gambar 7. Hasil simulasi model dua sel saraf terkopel dengan I = 0μA ........................ 12 Gambar 8. Hasil simulasi model dua sel saraf terkopel dengan I = 0,4μA ..................... 12 Gambar 9. Hasil simulasi model dua sel saraf terkopel dengan I = 0,41μA ................... 12 Gambar 10. Hasil simulasi model tiga sel saraf terkopel dengan I = 0μA........................ 12 Gambar 11. Hasil simulasi model tiga sel saraf terkopel dengan I = 0,4μA..................... 13 Gambar 12. Hasil simulasi model tiga sel saraf terkopel dengan I = 0,47μA................... 13 Gambar 13. Hasil simulasi model empat sel saraf terkopel dengan I = 0μA .................... 13 Gambar 14. Hasil simulasi model empat sel saraf terkopel dengan I = 0,5μA ................. 13 Gambar 15. Hasil simulasi model empat sel saraf terkopel dengan I = 0,54μA ............... 13 Gambar 16. Hasil simulasi model delapan sel saraf terkopel dengan I = 0μA ................. 13 Gambar 17. Hasil simulasi model delapan sel saraf terkopel dengan I = 0,55μA ............ 14 Gambar 18. Hasil simulasi model delapan sel saraf terkopel dengan I = 0,81μA ............ 14 Gambar 19. Hasil simulasi model enam belas sel saraf terkopel dengan I = 0μA ............ 14 Gambar 20. Hasil simulasi model enam belas sel saraf terkopel dengan I = 1μA ............ 14 Gambar 21. Hasil simulasi model enam belas sel saraf terkopel dengan I = 1,4μA ......... 14
DAFTAR LAMPIRAN halaman Lampiran 1. Diagram Alir Penelitian .................................................................
17
Lampiran 2. Rencana Kegiatan Penelitian .........................................................
18
Lampiran 3. Rencana Anggaran Penelitian ........................................................
18
Lampiran 4. Source Code Program dengan software Matlab .............................
18
BAB 1 PENDAHULUAN 1.1
Latar Belakang Sistem saraf merupakan salah satu sistem koordinasi yang berfungsi untuk menyampaikan rangsangan yang dideteksi dan direspon oleh tubuh. Sistem saraf memungkinkan makhluk hidup tanggap dengan cepat terhadap perubahan-perubahan yang terjadi di lingkungan luar maupun dalam. Sistem saraf terdiri atas banyak sel saraf, yang biasa disebut dengan neuron. Sel saraf berfungsi mengirimkan pesan yang berupa rangsangan atau tanggapan dari badan sel menuju ke dendrit. Setelah menerima berjuta-juta rangsangan informasi yang berasal baik dari luar maupun dari dalam tubuh, rangsangan tersebut diintegrasikan dan kemudian digunakan untuk menentukan respon apa yang akan diberikan oleh tubuh.1 Penjalaran atau propagasi serta proses integrasi impuls sel saraf merupakan hal yang menarik untuk dipelajari. Model matematika dinamika impuls pada satu sel saraf pertama kali dinyatakan oleh model Hodgkin-Huxley dan model FitzHugh-Nagumo.2,3,4 Model HogkinHuxley berbentuk sistem persamaan diferensial biasa nonlinear dengan empat variabel yang menggambarkan dinamika dari sel saraf, sedangkan model FitzHughNagumo mempunyai bentuk yang lebih sederhana, yaitu dalam bentuk sistem persamaan diferensial biasa dengan dua variabel yang otonom (autonomous).5,6 Pada tugas akhir ini akan dibuat sebuah model matematika untuk sel saraf terkopel dengan cara memodifikasi model sel saraf dari FitzHugh-Nagumo. Untuk mengetahui model sel saraf terkopel yang dibuat tersebut memiliki perilaku yang mendekati perilaku sel saraf sebenarnya, maka terlebih dahulu dibuat simulasi model satu sel saraf FitzHugh-Nagumo. Selanjutnya dibuat simulasi model sel saraf terkopel untuk mendapatkan hasil yang lebih baik. Hasil simulasi satu sel saraf FitzHughNagumo dengan sel saraf terkopel tersebut dibandingkan. Apabila perilakunya sudah sesuai, maka dilanjutkan dengan menganalisis arti fisis dari model sel saraf terkopel tersebut.
1.2 Tujuan Penelitian Penelitian ini bertujuan untuk: 1. Membuat model matematika untuk sel saraf terkopel. 2. Membuat simulasi dari model sel saraf terkopel yang telah dibuat. 3. Mempelajari dan menganalisis penjalaran atau propagasi impuls pada model sel saraf terkopel yang telah dibuat. 1.3 Perumusan Masalah 1. Bagaimana hasil simulasi dari model sel saraf terkopel yang telah dibuat ? 2. Apakah manfaat dari model tersebut ? 1.4 Hipotesis 1. Hasil simulasi dari model sel saraf terkopel lebih mendekati perilaku sel saraf sebenarnya dibandingkan dengan hasil simulasi satu sel saraf. 2. Model ini dapat digunakan untuk memperkirakan arus eksternal minimal yang dapat direspon oleh sel saraf terkopel.
BAB 2 TINJAUAN PUSTAKA 2.1
Sistem Saraf Sistem saraf tersusun oleh berjutajuta sel saraf yang mempunyai bentuk bervariasi. Sistem ini meliputi sistem saraf pusat dan sistem saraf tepi. Dalam kegiatannya, saraf mempunyai hubungan kerja seperti mata rantai (berurutan) antara receptor dan efector. Receptor adalah satu atau sekelompok sel saraf dan sel lainnya yang berfungsi mengenali rangsangan tertentu yang berasal dari luar atau dari dalam tubuh. Efector adalah sel atau organ yang menghasilkan tanggapan terhadap rangsangan, sebagai contoh: otot dan kelenjar.5 Sistem saraf terdiri dari jutaan sel saraf. Fungsi sel saraf adalah mengirimkan pesan (impuls) yang berupa rangsang atau tanggapan. Setiap sel saraf terdiri dari satu badan sel yang di dalamnya terdapat sitoplasma dan inti sel seperti pada Gambar 1. Dari badan sel keluar dua macam serabut saraf, yaitu dendrit dan akson (neurit). Setiap sel saraf hanya mempunyai satu akson dan minimal satu dendrit. Kedua serabut saraf ini berisi plasma sel.
2
Gambar 1. Bagian sel saraf 7
Pada bagian luar akson terdapat lapisan lemak disebut myelin yang merupakan kumpulan sel schwann yang menempel pada akson. Sel schwann adalah sel yang membentuk selubung lemak di seluruh serabut saraf myelin. Membran plasma sel schwann disebut neurilemma. Fungsi myelin adalah melindungi dan memberi nutrisi pada akson. Bagian dari akson yang tidak terbungkus mielin disebut nodus ranvier.1 Model Hodgkin-Huxley Alan Lloyd Hodgkin dan Andrew Fielding Huxley 2,8,9 melakukan percobaan yang bertujuan untuk menghitung fluks ionik dan perubahan permeabilitas membran pada mekanisme molekuler. Kemudian pada percobaan tersebut dikembangkan sebuah model deskripsi kinetik empiris yang cukup sederhana untuk membuat perhitungan praktis dari respon elektrik, namun cukup sesuai untuk memprediksi konduksi dari suatu rangsangan. Model tersebut tidak hanya terdiri dari persamaan matematika tetapi juga menunjukkan fitur utama dari mekanisme gerbang membran. Model tersebut kemudian disebut sebagai model Hodgkin-Huxley (model HH).
yang hanya dapat dimasuki oleh ion tertentu. Dua saluran utama membran yaitu saluran potasium (ion K+) dan saluran sodium (ion Na+) yang bisa menyebabkan potensial aksi, seperti yang terlihat pada Gambar 2. Kedua saluran tersebut secara hipotesis dipengaruhi oleh tiga variabel yaitu m, h dan n. Variabel m merupakan peluang terbukanya satu rapidly activating channels untuk natrium, sedangkan h adalah peluang terbukanya satu slowly activating channels untuk kalium dan variabel n yang merupakan peluang terbukanya satu slowly activating channels untuk natrium.
2.2
Gambar 2. Skema potensial aksi 10
Dari grafik hasil simulasi parameter Hodgkin-Huxley 11 diketahui bahwa m berubah secara cepat sehingga berpengaruh terhadap terjadinya depolarisasi (depolarization) potensial aksi. Kemudian diikuti dengan variabel n yang menyebabkan potensial aksi mencapai puncak, Pada saat ini, potensial aksi mengalami hyperpolarization. Variabel h merupakan variabel yang paling lambat berubah, Variable ini memiliki peranan sebagai penyebab potensial aksi turun melebihi batas potensial istirahat. Sistem persamaan empat variabel dari Hodgkin-Huxley 2,8,12 adalah sebagai berikut:
Model HH menjelaskan bahwa pada membran terdapat saluran-saluran khusus
dV 1 m 3 hg Na E E Na n 4 g K E E K g L E E L I ................................ (1) dt 3 dn n 1 n n n ............................................................................................................ (2) dt
dm m 1 m m m ........................................................................................................ (3) dt dh h 1 h h h ............................................................................................................ (4) dt
3
Keterangan : V adalah potensial membran. I adalah arus ionik total yang melewati membran. m adalah probabilitas salah satu dari tiga partikel aktivasi yang diperlukan untuk berkontribusi dalam aktivasi gerbang Na. m3 adalah probabilitas ketiga partikel aktivasi yang telah menghasilkan sebuah channel terbuka. h adalah probabilitas satu partikel yang tidak aktif yang tidak menyebabkan gerbang Na menutup. adalah konduktansi sodium maksimal. E adalah potensial membran total. ENa adalah potensial membran Na. n adalah probabilitas salah satu dari empat partikel aktivasi yang telah mempengaruhi keadaan dari gerbang K. adalah konduktansi potassium maksimal. EK adalah potensial membran K. adalah maksimal leakage conductance. EL adalah potensial membran leakage. αm adalah kontanta untuk partikel yang tidak mengaktifkan gerbang. βm adalah konstanta untuk partikel yang mengaktifkan gerbang. 2.3
Model FitzHugh-Nagumo Richard FitzHugh dan Nagumo membuat suatu model dengan cara menyederhanakan sistem empat variabel model HH tersebut menjadi sistem dua variabel sehingga lebih mudah untuk dianalisa. Sifat model ini mirip dengan model Hodgkin-Huxley secara kualitatif. Model FitzHugh-Nagumo dapat diaproksimasikan secara kualitatif. Bentuk persamaannya adalah sebagai berikut:
dx f x y I ........................(5) dt dy bx y .................................(6) dt
f x xa x x 1 ..........................(7)
dengan 0 < a <1 dan b, adalah konstanta positif.5 Berdasarkan aproksimasi ini FitzHugh 3 membuat persamaan yang berasal dari persamaan umum osilasi teredam, yaitu: d 2x dx k x 0 ............................. (8) 2 dt dt k merupakan konstanta redaman
FitzHugh menggunakan persamaan Van der Pool 13 yaitu merubah konstanta damping dengan fungsi kuadrat:
d 2x dx c x2 1 x 0 .......... (9) 2 dt dt c merupakan konstanta positif. Untuk mendapatkan sistem dua variabel dilakukan transformasi dengan metode Lienard 14, dan didapatkan model FitzHugNagumo sebagai berikut : ........ (10) dx x3 c y x I dt 3
dy x a by c ............. (11) dt Keterangan : a dan b merupakan konstanta positif. Variabel x berperan sebagai potensial membran V pada model HH Variabel y berperan sebagai variabel m, n, dan h pada model HH 2.4
Model Matematika Sel Saraf Terkopel Model saraf terkopel yang penulis gunakan ini berdasarkan model matematika dari Hindmarsh-Rose 15, dimana pada model Hindmarsh-Rose menambahkan fungsi kopel pada salah satu sistem persamaan diferensialnya, seperti terlihat di bawah ini:
n
h xi , x j xi Vs g s cij x j j 1
.................... (12) Dimana merupakan kopel dari sinapsis yang dimodelkan dengan fungsi sigmoid dengan memiliki ambang batas potensial aksi di setiap sel saraf.
x j
1
1 exp x j s
........ (13)
adalah batas ambang yang dicapai setiap potensial aksi sel saraf. Sel saraf seharusnya identik dan sinapsis yang cepat. Pada model Hindmarsh-Rose parameter = 0.17 sesuai dengan kekuatan kopling dari sinapsis. Karena sinapsis merupakan exitatory maka potensial pembalik ( ) harus lebih besar daripada untuk semua nilai i dan t. Dari pengetahuan ini, model matematika yang akan dibuat adalah sebagai berikut:
4
n dxi x3 c yi xi i I xi Vs g s cij x j ..................................... (14) dt 3 j 1 dy i xi a by i c .................................................................................. (15) dt
x j
1 ........................................................................................ (16) 1 exp s x j
adalah matrik penghubung yang elemennya memenuhi syarat sebagai berikut: jika i dan j saling terhubung jika i dan j tidak saling terhubung i = 1,2,3...,n, j = 1,2,3....,n Pada model ini diasumsikan arah perambatan impuls itu bidirection yaitu impuls bisa merambat dua arah. Misalkan rangsangan terjadi pada bagian tengah kumpulan sel saraf, maka arah impuls bisa ke kanan dan ke kiri. Tetapi fokus perambatan yang dipelajari dalam model ini hanya dari efector menuju receptor. Seperti perambatan dari ujung garis (dalam hal ini efector) yang tersusun dari kumpulan beberapa sel saraf ke ujung satunya (receptor ). 2.5 Mekanisme Penghantar Impuls 2.5.1 Penghantaran Impuls Melalui Sel Saraf Penghantaran impuls yang berupa rangsangan atau tanggapan melalui serabut saraf dapat terjadi karena adanya perbedaan potensial listrik antara bagian luar dan bagian dalam sel saraf. Pada waktu sel saraf beristirahat, kutub positif terdapat di bagian luar dan kutub negatif terdapat di bagian dalam sel saraf. Rangsangan pada indra manusia menyebabkan terjadinya pembalikan perbedaan potensial listrik sesaat, perubahan potensial ini terjadi berurutan sepanjang serabut saraf. Bila impuls telah lewat maka untuk sementara serabut sel saraf tidak dapat dilalui oleh impuls lagi, karena terjadi perubahan potensial kembali seperti semula (potensial istirahat). Untuk dapat berfungsi kembali diperlukan waktu sekitar 0,001 sampai 0,002 detik.16 Rangsangan yang kurang kuat atau dibawah batas ambang (threshold) tidak akan menghasilkan impuls yang dapat merubah potensial listrik, atau dengan kata lain tubuh tidak merespon rangsangan tersebut. Tetapi bila kekuatan rangsangan
tersebut di atas threshold maka impuls akan dihantarkan sampai ke ujung sel saraf. 2.5.2 Penghantaran Impuls Melalui Sinapsis Titik pertemuan antara terminal akson salah satu sel saraf dengan sel saraf lain dinamakan sinapsis. Setiap terminal akson membengkak membentuk tonjolan sinapsis. Di dalam sitoplasma tonjolan sinapsis terdapat struktur kumpulan membran kecil berisi neurotransmitter yang disebut vesikula sinapsis. Sel saraf yang berakhir pada tonjolan sinapsis disebut neuron pra-sinapsis. Membran ujung dendrit dari sel berikutnya yang membentuk sinapsis disebut post-sinapsis. Bila impuls sampai pada ujung sel saraf, maka vesikula bergerak dan melebur dengan membran prasinapsis. Kemudian vesikula akan melepaskan neurotransmitter berupa asetilkolin. Neurontransmitter adalah suatu zat kimia yang dapat menyeberangkan impuls dari neuron pra-sinapsis ke postsinapsis. Neurontransmitter ada bermacammacam misalnya asetilkolin yang terdapat di seluruh tubuh, noradrenalin terdapat di sistem saraf simpatik, dan dopamin serta serotonin yang terdapat di otak. Asetilkolin berdifusi melewati celah sinapsis dan menempel pada receptor yang terdapat pada membran post- sinapsis. Penempelan asetilkolin pada receptor menimbulkan impuls pada sel saraf berikutnya. Bila asetilkolin sudah melaksanakan tugasnya maka akan diuraikan oleh enzim asetilkolinesterase yang dihasilkan oleh membran post-sinapsis.16 2.6
Metode Ode45 Persamaan diferensial biasa pada MATLAB bisa diselesaikan dengan metode ode45. Metode ini merupakan implementasi dari Runge-Kutta dengan variabel waktu sebagai iterasi. Metode Runge-Kutta 17 secara umum sebagai berikut :
........................................ (17) Pendekatan solusi dengan menggunakan metode Runge-Kutta orde-4 adalah : .......................................................... (18) Sedangkan untuk mendapatkan nilai yang lebih baik untuk solusi tersebut digunakan metode Runge-Kutta orde 5 sebagai berikut: ......................................................... (19)
[t, y] = ode45(„fname‟, tspan, y0); Keterangan : Fname : Nama fungsi dari Mfile yang digunakan, biasanya didefinisikan sebagai berikut : function dydt = fname(t, y) tspan : dua elemen vektor yang mendefinisikan rentang dari waktu awal dan waktu akhir. y0 : vektor dari kondisi awal untuk variabel y
[t, s] = ode45(„subinagumo‟, [1 30], [0 0]); Plot (t, s(:,1)); Hasil simulasinya adalah sebagai berikut : 2 1.5 1 Potensial aksi(mV)
Sedangkan untuk metode ode45 18 sintak umum pada MATLAB sebagai berikut :
Contoh program model satu sel saraf sebagai berikut : function dsdt =subinagumo(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I = -0.35; x = s(1); % kondisi awal x y = s(2); % kondisi awal y % persamaan diferensialnya dsdt(1) = c*(x+y-x^3/3+I); % pers 22 dsdt(2) = -(x-a+b*y)/c; % pers 23 Simpan dengan nama file subinagumo.m, untuk mejalankan program tersebut ketik pada command window :
0.5 0 -0.5 -1 -1.5 -2
0
5
10
15 Waktu(ms)
20
25
30
BAB 3 METODOLOGI PENELITIAN 3.1
Waktu dan Tempat Penelitian Penelitian ini dilakukan di Laboratorium Fisika Teori dan Komputasi, Departemen Fisika, FMIPA, IPB dari bulan September 2010 sampai Januari 2011.
3.2
3.3
Metode Penelitian Metode penelitian ini adalah pembuatan sebuah program simulasi sederhana dari model satu sel saraf FitzHugh-Nagumo (persamaan 10 dan 11) menggunakan software MATLAB R2010b. Selanjutnya membuat model matematika sel saraf terkopel beserta simulasinya. Setelah itu simulasi ini dibandingkan dengan simulasi satu sel saraf FitzHugh-Nagumo, apabila perilakunya sudah sama selanjutnya dianalisis sifat fisis dari penjalaran impuls pada sel saraf terkopel tersebut.
BAB 4 HASIL DAN PEMBAHASAN 4.1
Hasil Simulasi Model Satu Sel Saraf Hasil simulasi model satu sel saraf dari Fitzhugh-Nagumo untuk I = 0 samapai 0,35 μA dapat dilihat pada Gambar 3 sampai 7. Nilai parameter a = 0,7; b = 0,8; dan c = 3, nilai parameter ini dipilih sama dengan nilai pada simulasi yang dilakukan Fitzhugh.3 1.8 1.6 1.4
Potensial aksi(mV)
Peralatan Peralatan yang digunakan dalam penelitian ini adalah sebuah laptop dengan processor AMD Turion X-2 (2.2 GHz), 3 GB RAM. Software yang digunakan dalam penelitian ini adalah MS. Office 2007 dan MATLAB R2010b. Pada penelitian ini penulis menggunakan sumber pustaka berupa jurnal-jurnal ilmiah neuron dan buku-buku tentang neuron.
1.2 1 0.8 0.6 0.4 0.2
3.3.2 Analisa numerik Analisa ini dilakukan karena model yang didapatkan pada perancangan model sel saraf ini merupakan sistem dua persamaan diferensial autonomous, dimana pada persamaan ini sangat sulit untuk diselesaikan secara analitik. Model yang dibuat berbentuk persamaan diferensial biasa, sehingga metode numerik yang akurat adalah Runge Kutta orde 45 atau ode45. 3.3.3 Analisa sifat fisis model yang telah dibuat Analisa ini dilakukan dengan melihat hasil simulasi yang diperoleh dari model yang telah dibuat tersebut. Ini dilakukan untuk mengetahui penjalaran impuls yang sebenarnya. Selain itu juga digunakan untuk memperkirakan arus eksternal yang dapat direspon oleh tubuh.
0
0
5
10
15 Waktu(ms)
20
25
30
Gambar 3. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0 μA.
Gambar 3 menunjukan simulasi sebelum ada arus eskternal, terlihat tidak ada impuls yang menjalar sehingga berakibat potensial aksi dari sel saraf lama-kelamaan akan menuju ke keadaan istirahat atau sering di sebut potensial istirahat. Hal ini dapat dijelaskan seperti tubuh manusia yang tidak mendapatkan rangsang dari luar maka tubuh juga tidak akan melakukan respon apapun. Oleh karena itu, simulasi ini sudah bisa menjelaskan keadaan normal dari tubuh manusia yang belum mendapatkan rangsang dari luar. 2
1.5
1 Potensial aksi(mV)
3.3.1 Studi pustaka Studi pustaka dilakukan untuk memahami proses dinamika impuls dalam sel saraf sehingga memudahkan perancangan program simulasinya. Kemudian melihat hubungan antara grafik-grafik yang akan dihasilkan dalam simulasi dengan sifat fisis dari sel saraf itu sendiri. Adanya studi pustaka ini akan membantu penulis dalam menganalisa hasil yang di dapat dalam simulsi sel saraf tersebut.
0.5
0
-0.5
-1
-1.5 0
5
10
15 Waktu(ms)
20
25
Gambar 4. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,2 μA.
30
7
2 1.5
Potensial aksi(mV)
1 0.5 0 -0.5 -1 -1.5 -2
0
5
10
15 Waktu(ms)
20
25
30
Gambar 5. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,3 μA. 2 1.5
Potensial aksi(mV)
1 0.5 0 -0.5 -1 -1.5 -2
0
5
10
15 Waktu(ms)
20
25
Gambar 6. Hasil simulasi model satu sel saraf FitzHugh-Nagumo dengan I = 0,35 μA.
Ketika arus eksternal diterapkan ke dalam sel saraf tubuh, impuls (rangsangan) tersebut akan direspon oleh sel saraf tubuh atau tidak direspon oleh sel saraf tubuh, tergantung dari berapa besar arus yang diterapkan. Seperti seseorang yang dipukul sangat keras pasti akan memberikan respon
30
pada rangsangan itu, tetapi lain halnya jika orang yang tertepa angin yang sangat pelan belum tentu orang tersebut merespon rangsangan tersebut. Atau dengan kata lain ada batas minimal dari arus yang dapat direspon oleh seseorang. Berdasarkan hasil simulasi ini terlihat bahwa ketika arus yang diterapkan itu sebesar 0.2 μA maka penjalaran impuls itu belum terjadi seperti pada Gambar 4. Ketika arus yang diterapkan itu 0.3 μA sudah ada sedikit osilasi tapi masih dalam bentuk osilasi teredam seperti Gambar 5, artinya rangsangan ini belum mendapat respon dari tubuh. Tetapi ketika arusnya sebesar 0.35 μA mulai terjadi osilasi yang tidak teredam, hal ini menunjukan impuls itu mulai menjalar dan rangsangan sebesar ini dapat direspon oleh tubuh. Tetapi sebenarnya yang mempengaruhi respon terhadap rangsang itu tidak hanya arus dari luar tubuh saja, kondisi tubuh juga sangat berpengaruh, misalkan orang yang sakit dengan yang sehat akan mengalami respon yang berbeda, tetapi dalam simulasi ini hanya dibahas pengaruh arus saja. Mungkin untuk penelitian lebih lanjut bisa ditambahkan faktor-faktor tersebut. 4.2
Model Matematika Sel Saraf Terkopel 4.2.1 Model Dua Sel Saraf Terkopel Matriks
0 1 C2 1 0
Sistem persamaan diferensial terkopel yang terbentuk:
3 dx1 x1 1 c y x I 1 1 x1 Vs 0.17 1 exp x dt 3 s 2 dy1 x1 a by1 c dt 3 dx2 x 1 c y 2 x2 2 I x2 Vs 0.17 1 exp x dt 3 s 1 dy 2 x 2 a by 2 c ................................................................................................ (20) dt
8
4.2.2 Model Tiga Sel Saraf Terkopel
0 1 1 Matriks C 3 1 0 1 1 1 0 Sistem persamaan diferensial terkopel yang terbentuk : 3 dx1 x1 1 1 c y x I 1 1 x1 Vs 0.17 1 exp x 1 exp x dt 3 s 2 s 3
dy1 x1 a by1 c dt 3 dx2 x 1 1 c y 2 x2 2 I x2 Vs 0.17 dt 3 1 exp s x1 1 exp s x3
dy 2 x 2 a by 2 c dt 3 dx3 x 1 1 c y3 x3 3 I x3 Vs 0.17 dt 3 1 exp x 1 exp x s 1 s 2
dy 3 x3 a by 3 c ........................................................................................... (21) dt 4.2.3 Model Empat Sel Saraf Terkopel
0 1 Matriks C 4 1 1
1 0 1 1
1 1 0 1
1 1 1 0
Sistem persamaan diferensial terkopel yang terbentuk : 3 dx1 x 1 1 1 c y1 x1 1 I x1 Vs 0.17 dt 3 1 exp x 1 exp x 1 exp x s 2 s 3 s 4
dy1 x1 a by1 c dt 3 dx2 x 1 1 1 c y 2 x 2 2 I x 2 Vs 0.17 dt 3 1 exp s x1 1 exp s x3 1 exp s x 4
dy 2 x 2 a by 2 c dt
9
3 dx3 x 1 1 1 c y 3 x3 3 I x3 Vs 0.17 dt 3 1 exp x 1 exp x 1 exp x s 1 s 2 s 4
dy 3 x3 a by 3 c dt 3 dx4 x 1 1 1 c y 4 x4 4 I x4 Vs 0.17 dt 3 1 exp x 1 exp x 1 exp x s 1 s 2 s 3 dy 4 x 4 a by 4 c dt ........................................................................................... (22)
4.2.4 Model Delapan Sel Saraf Terkopel
0 1 1 1 Matriks C 8 1 1 1 1
1 0 1 1 1 1 1 1
1 1 0 1 1 1 1 1
1 1 1 0 1 1 1 1
1 1 1 1 0 1 1 1
1 1 1 1 1 0 1 1
1 1 1 1 1 1 0 1
1 1 1 1 1 1 1 0
Sistem persamaan diferensial terkopel yang terbentuk : 1 1 1 1 exp s x 2 1 exp s x3 1 exp s x 4 3 dx1 x 1 1 1 c y1 x1 1 I x1 Vs 0.17 dt 3 1 exp s x5 1 exp s x6 1 exp s x 7 1 1 exp x s 8
dy1 x1 a by1 c dt 1 1 1 1 exp s x1 1 exp s x3 1 exp s x 4 3 dx 2 x2 1 1 1 c y 2 x 2 I x 2 Vs 0.17 dt 3 1 exp s x5 1 exp s x 6 1 exp s x 7 1 1 exp x s 8
dy 2 x 2 a by 2 c dt
10
1 1 1 1 exp s x1 1 exp s x 2 1 exp s x 4 3 dx3 x3 1 1 1 c y 3 x 3 I x3 Vs 0.17 1 exp s x5 1 exp s x 6 1 exp s x 7 dt 3 1 1 exp x s 8s
dy 3 x3 a by 3 c dt
1 1 1 1 exp s x1 1 exp s x 2 1 exp s x3 3 dx8 x 1 1 1 c y8 x8 8 I x8 Vs 0.17 dt 3 1 exp s x 4 1 exp s x5 1 exp s x 6 1 1 exp x s 7 dy 8 x8 a by 8 c dt ............................................................................................ (23)
4.2.5
Model Enam Belas Sel Saraf Terkopel
Matriks C16 adalah sebagai berikut :
C16
0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
11
Sistem persamaan diferensial terkopel yang terbentuk : 1 1 exp s 1 1 exp s 3 dx1 x 1 c y1 x1 1 I x1 V s 0.17 1 exp s dt 3 1 1 exp s 1 1 exp s
1 1 1 exp s x 3 1 exp s x 4 1 1 x 5 1 exp s x 6 1 exp s x 7 1 1 x8 1 exp s x 9 1 exp s x10 1 1 x11 1 exp s x12 1 exp s x13 1 1 x14 1 exp s x15 1 exp s x16
x2
dy1 x1 a by1 c dt 1 1 exp s 1 1 exp s 3 dx 2 x 1 c y 2 x 2 2 I x 2 V s 0.17 dt 3 1 exp s 1 1 exp s 1 1 exp s
1 1 1 exp s x 3 1 exp s x 4 1 1 x 5 1 exp s x 6 1 exp s x 7 1 1 x8 1 exp s x 9 1 exp s x10 1 1 x11 1 exp s x12 1 exp s x13 1 1 x14 1 exp s x15 1 exp s x16 x1
dy 2 x 2 a by 2 c dt 1 1 exp s 1 1 exp s 3 dx 3 x 1 c y 3 x 3 3 I x 3 V s 0.17 dt 3 1 exp s 1 1 exp s 1 1 exp s
1 1 x 5 1 exp s x 6 1 exp s x 7 1 1 x8 1 exp s x 9 1 exp s x10 1 1 x11 1 exp s x12 1 exp s x13 1 1 x14 1 exp s x15 1 exp s x16
x1
1 1 1 exp s x 2 1 exp s x 4
dy 3 x3 a by 3 c dt
1 1 1 1 exp s x1 1 exp s x 2 1 exp s x3 1 1 1 3 dx16 x 1 exp s x 4 1 exp s x8 1 exp s x9 c y16 x16 16 I x16 Vs 0.17 1 1 1 dt 3 1 exp s x10 1 exp s x11 1 exp s x12 1 1 1 1 exp x 1 exp x 1 exp x s 13 s 14 s 15 dy16 x16 a by16 c dt .................................................................................... (24)
12
4.3
1.8 1.6
1.2 1 0.8 0.6 0.4 0.2 0
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 7. Hasil simulasi model dua sel saraf terkopel dengan I = 0 μA. 2 1.5
Potensial aksi(mV)
1
1.5 1 0.5 0 -0.5 -1 -1.5 -2
0
5
10
20
25 30 Waktu(ms)
35
40
45
50
Hasil simulasi dua sel saraf terkopel yang telah dibuat memiliki bentuk grafik yang sama dengan model satu sel saraf FitzHugh-Nagumo. Hal ini menunjukkan bahwa model yang dibuat sudah bisa memenuhi perilaku sel saraf. Oleh karena itu, selanjutnya dilakukan analisa hasil simulasinya. Hasil simulasi model dua sel saraf terkopel ini dapat dilihat pada Gambar 7 sampai 9, dimana ketika tidak ada arus eksternal impuls tidak menjalar seperti pada Gambar 7. Ketika arus yang diterapkan 0,4 μA, impuls sudah terjadi sedikit osilasi teredam. Artinya rangsangan dengan arus sebesar ini belum mendapatkan respon dari tubuh. Ketika arus yang diterapkan 0,41 μA impuls mulai menjalar dan hal ini menunjukan bahwa rangsangan dengan arus sebesar itu telah direspon oleh tubuh. 4.3.2
0.5
15
Gambar 9. Hasil simulasi model dua sel saraf terkopel dengan I = 0,41 μA.
Hasil Simulasi Tiga Sel Saraf Terkopel
0
1.8
-0.5
1.6 -1
1.4 -1.5 -2
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 8. Hasil simulasi model dua sel saraf terkopel dengan I = 0,4 μA.
Potensial aksi(mV)
Potensial aksi(mV)
1.4
2
Potensial aksi(mV)
Hasil Simulasi Model Sel Saraf Terkopel Hasil simulasi model sel saraf terkopel ( Persamaan 20 sampai 24 ) dapat dilihat pada Gambar 7 sampai 9 untuk model dua sel saraf terkopel, Gambar 10 sampai 12 untuk model tiga sel saraf terkopel, Gambar 13 sampai 15 untuk model empat sel saraf terkopel, Gambar 16 sampai 18 untuk model delapan sel saraf terkopel, Gambar 19 sampai 21 untuk model enam belas sel saraf terkopel. Nilai parameter a=0,7; b=0,8; c = 3; =2 mV; dan =-0,25 mV, nilai parameter ini diambil berdasarkan simulasi yang dilakukan Fitzhugh 3 dan Belykh.17 4.3.1 Hasil simulasi dua sel saraf terkopel
1.2 1 0.8 0.6 0.4 0.2 0
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 10. Hasil simulasi model tiga sel saraf terkopel dengan I = 0 μA.
13
2
2
1.5
1.5 1
0.5
Potensial aksi(mV)
Potensial aksi(mV)
1
0 -0.5
0.5 0 -0.5
-1
-1
-1.5 -2
-1.5
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
-2
50
Gambar 11. Hasil simulasi model tiga sel saraf terkopel dengan I = 0,4 μA.
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 14. Hasil simulasi model empat sel saraf terkopel dengan I = 0,5 μA
2 2 1.5 1.5 1
Potensial aksi(mV)
Potensial aksi(mV)
1 0.5 0 -0.5
0.5 0 -0.5
-1 -1 -1.5 -1.5 -2
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
-2
Gambar 12. Hasil simulasi model tiga sel saraf terkopel dengan I = 0,47 μA.
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 15. Hasil simulasi model empat sel saraf terkopel dengan I = 0,54 μA.
Gambar 10 menunjukan tidak ada arus eksternal dan impuls belum menjalar. Ketika arus yang diterapkan 0,4 μA, impuls belum menjalar. Artinya rangsangan dengan arus sebesar ini belum mendapatkan respon dari tubuh.Ketika arus yang diterapkan 0,47 μA impuls mulai menjalar seperti Gambar 12. Hal ini menunjukan bahwa rangsangan dengan arus sebesar itu telah direspon oleh tubuh. 4.3.3 Hasil Simulasi Empat Sel Saraf Terkopel
Gambar 13 menunjukan tidak ada arus eksternal dan impuls belum menjalar. Ketika arus yang diterapkan 0,5 μA, impuls belum menjalar.Ketika arus yang diterapkan 0,54 μA impuls mulai menjalar seperti Gambar 15. Hal ini menunjukan bahwa rangsangan dengan arus sebesar itu telah direspon oleh tubuh. 4.3.4 Hasil Simulasi Delapan Sel Saraf Terkopel 1.8
1.8
1.6
1.6
1.4
Potensial aksi(mV)
Potensial aksi(mV)
1.4 1.2 1 0.8
1.2 1 0.8
0.6
0.6
0.4
0.4
0.2 0
0.2 0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 13. Hasil simulasi model empat sel saraf terkopel dengan I = 0 μA.
0
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 16. Hasil simulasi model delapan sel saraf terkopel dengan I = 0 μA.
14
1.4
2
1.2
1.5
1 Potensial aksi(mV)
Potensial aksi(mV)
1
0.5
0
0.8 0.6 0.4
-0.5
0.2
-1
0
-1.5 0
-0.2 5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 17. Hasil simulasi model delapan sel saraf terkopel dengan I = 0,55 μA.
20
25 30 Waktu(ms)
35
40
45
50
1
0.5
Potensial aksi(mV)
Potensial aksi(mV)
15
1.5
1
0 -0.5 -1
0.5 0 -0.5 -1
-1.5
-1.5 0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 18. Hasil simulasi model delapan sel saraf terkopel dengan I = 0,81 μA.
Gambar 16 menunjukan tidak ada arus eksternal dan impuls belum menjalar. Ketika arus yang diterapkan 0,55μA, impuls belum menjalar. Ketika arus yang diterapkan 0,81 μA impuls mulai menjalar seperti Gambar 18. Hal ini menunjukan bahwa arus minimal yang dapat direspon oleh tubuh adalah 0,81 μA. 4.3.5 Hasil Simulasi Enam Belas Sel Saraf Terkopel 1.8 1.6 1.4
Potensial aksi(mV)
10
2
1.5
1.2 1 0.8 0.6 0.4 0.2 0
5
Gambar 20. Hasil simulasi model enam belas sel saraf terkopel dengan I = 1 μA.
2
-2
0
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
50
Gambar 19. Hasil simulasi model enam belas sel saraf terkopel dengan I = 0 μA.
-2
0
5
10
15
20
25 30 Waktu(ms)
35
40
45
Gambar 21. Hasil simulasi model enam belas sel saraf terkopel dengan I = 1,4 μA.
Gambar 19 menunjukan tidak ada arus eksternal dan impuls belum menjalar. Ketika arus yang diterapkan 1 μA, impuls belum menjalar.Ketika arus yang diterapkan 1,4 μA impuls mulai menjalar seperti Gambar 21. Hal ini menunjukan bahwa rangsangan dengan arus sebesar itu telah direspon oleh tubuh. 4.4 Analisa Model Sel Saraf Terkopel Berdasarkan simulasi sel saraf terkopel yang diperoleh terlihat bahwa semakin banyak sel saraf terkopel maka arus minimal yang dapat direspon tubuh juga semakin bertambah besar. Dari model yang telah dibuat didapatkan hasil simulasi untuk menentukan arus minimal yang dapat direspon oleh tubuh. Dari hasil simulasi dua sel saraf terkopel sampai enam belas sel saraf terkopel terlihat kenaikan arus minimal yang dapat direspon tubuh itu sekitar 0,06 μA. Hal ini bisa digunakan untuk memperkirakan arus minimal yang dapat direspon tubuh pada sel saraf sebenarnya. Misalkan pada tubuh jumlah sel saraf yang bekerja pada penghataran impuls ada satu juta, maka dapat dicari arus minimalnya agar dapat direspon oleh tubuh. Hal ini dapat
50
dianalogikan jika untuk satu sel saraf arus minimalnya 0,06 μA, maka untuk satu juta sel saraf sekitar 0,06 A. Jadi rangsangan yang datang dari luar tubuh yang memiliki arus sebesar 0,06 A akan direspon oleh tubuh, sedangkan rangsangan yang arusnya kurang dari itu tidak akan direspon oleh tubuh.
BAB 5 KESIMPULAN DAN SARAN 5.1
Kesimpulan Model matematika untuk sel saraf terkopel yang dibuat telah sesuai dengan perilaku sel saraf. Hal ini dapat dilihat dari hasil simulasi model sel saraf terkopel tersebut memiliki hasil simulasi yang sama dengan model satu sel saraf FitzHughNagumo. Model FitzHugh-Nagumo ini yang menjadi acuan model yang dibuat pada skripsi ini, karena model FitzHugh-Nagumo ini telah diakui oleh banyak ilmuwan. Simulasi satu sel saraf ini hanya bisa menjelaskan perilaku satu sel saraf, karena di dalam tubuh yang berperan dalam penjalaran impuls itu tidak hanya satu sel saraf saja, melainkan banyak sel saraf. Oleh karena itu, model sel saraf terkopel yang dibuat lebih bisa menjelaskan perilaku sel saraf sebenarnya dalam merespon rangsangan. Hasil simulasi satu sel saraf terkopel menunjukan bahwa arus minimal agar impuls dapat direspon yaitu 0,35 μA. Sedangkan hasil simulasi untuk sel saraf terkopel menujukan hasil arus yang semakin besar dengan semakin banyak model sel saraf yang digunakan. Untuk model dua sel saraf terkopel arus minimal agar impuls dapat direspon yaitu 0,41 μA, untuk model tiga sel saraf terkopel arus minimalnya yaitu 0,47 μA, untuk model empat sel saraf terkopel arus minimalnya yaitu 0,54 μA, untuk model delapan sel saraf terkopel membutuhkan arus minimal 0,81 μA agar bisa merespon impuls tersebut, sedangkan untuk model enam belas sel saraf terkopel membutuhkan arus minimal 1,4 μA untuk meresponnya. Hasil simulasi sel saraf terkopel ini menunjukan peningkatan arus minimal agar bisa merespon impuls tersebut. Dimana peningkatan arus ini sekitar 0,06 μA setiap
penambahan satu sel saraf. Oleh karena itu dengan menganggap bahwa jumlah sel saraf dalam tubuh yang berperan dalam penghantaran impuls itu sekitar satu juta sel saraf maka dapat diprediksi arus minimal agar impuls dapat direspon oleh tubuh sekitar 0,06 A, sedangkan arus yang kurang dari itu tidak akan mendapat respon dari tubuh. 5.2
Saran Simulasi ini hanya membahas mekanisme respon saraf secara umum, yaitu dalam merespon rangsang dari luar sel saraf dalam kondisi sama untuk semua tubuh manusia. Pada simulasi ini tidak memperhitungkan kondisi tubuh manusia itu, misalkan tubuh dalam keadaan sehat atau tidak. Karena kondisi tubuh dalam keadaan sehat dan kurang sehat akan merespon rangsang secara berbeda. Untuk penelitian selanjutnya bisa membuat simulasi model saraf dengan memperhatikan faktor-faktor tersebut. Sehingga akan didapatkan sebuah simulasi yang lebih komplek dan banyak mewakili peristiwa yang sering terjadi dalam kehidupan seharihari. Untuk melakukan simulasi ini minimal digunakan komputer dengan processor pentium 4, 1GB RAM. Untuk mendapatkan grafik yang lebih bagus pada MATLAB R2010b disarankan menggunakan komputer dengan processor diatas pentium 4, 2GB RAM. Hal ini selain bisa mendapatkan grafik yang lebih bagus juga memepercepat running program.
DAFTAR PUSTAKA 1.
2.
3.
4.
Oswari, S. (2008). Model Matematika Penjalaran Impuls Saraf pada Satu Sel Saraf di Subthalamik Nukleus. Bandung, 1-3. Hodgkin, A. L. & Huxley, A. F. (1952). A Quantitative Description Of Membranecurrent and Application to Conduction and Excitation in Nerve. J. Physiol, 117, 500-544. FitzHugh, R. (1961). Impuls and Physiology States in Theoritical Models of Nerve Membran. Biophysical, 1, 445-447. Izhikevich, E. M. (2002). Dynamical System in Neuroscience : The Geometry of Excitability and
16
5.
6. 7.
8.
9.
10.
11.
Brusting. London: The MIT Press Cambridge, 89-157. Yuanita, D. (2009). Dinamika Impuls pada Satu Sel Saraf dengan Akupuntur sebagai Stimulus. Bandung, 12-13. Alatas, H. (2009). Fisika Nonlinear : Dinamika Sistem. edisi 1. Bogor. Faisal. “ Jaringan Saraf Tiruan.” Blogspot. 09 Oktober 2010. Web. 20 Oktober 2010.
. Nagumo, J. (1962). An Active Pulse Transmission Line Simulating Nerve Axon. Proc. IRE, 501, 2061-2063. Hodgkin, A. L. & Huxley, A. F. (1952). The Components of Membrane Conductance in The Giant Axon of Loligo. J. Physiol, 116, 473496. Bay02pisay. “Nerve impulse transmission”. Wordpress. 15 Januari 2009. Web. 20 Oktober 2010. . Hodgkin, A. L. & Huxley, A. F. (1952). Currents Caried by Sodium
12.
13. 14.
15.
16.
17.
18.
and Potasium Ion Through The Membrane of The Giant Axon of Loligo. J. Physiol, 116, 449-472. Corson, N. & Aziz, A. (2006). Dynamics and Complexity of Hindmarsh-Rose Neuronal System. Encyclopedia of Mathematical Physics: Elsevier, 5, 213-226. Pol, V. D. (1926). On Relaxation Oscillation. Phil. Mag, 2, 978-979. Minorsky, N. (1947). Introduction to Nonlinear Mechanics, Edward Brother. Inc, 105-110. Lange, E., Belykh, I. & Hasler, M. (2004). Synchronization of Bursting Neurons: What matters in the Network Topology. EPLF, Switzerland. Anonim. “Sistem Saraf.” Web. 20 oktober 2010. Mathews, J. H. & Fink, K. D. (1999). Numerical Method Using Matlab. Prentice Hall, 458-474. Shampine, L. F., Gladwell, I. & Thompson, S. (2003). Solving ODEs with MATLAB. Cambridge University Press, NewYork.
LAMPIRAN Lampiran 1. Diagram Alir Penelitian
Mulai
Penelusuran literatur
Sudah siap
Pembuatan program simulasi satu sel saraf
Perancangan dan pembuatan model matematika untuk sel saraf terkopel
Aplikasi Matlab
Pembuatan program simulasi sel saraf terkopel
Hasil simluasi satu sel saraf
Aplikasi Matlab
Hasil simulasi sel saraf terkopel
Hasilnya dibandingkan
Tidak
Sesuai
Analisis
Penyusunan laporan
Selesai
18
Lampiran 2. Rencana Kegiatan Penelitian Bulan september 2010-Januari 2011
Kegiatan
Sep
Okt
Nov
Des
Jan
Penelusuran litertur Penyusunan proposal usulan penelitian dan kolokium Pembuatan program simulasi satu sel saraf Pembuatan program simulasi sel saraf terkopel Analisa akhir Penyusunan laporan
Lampiran 3. Rencana Anggaran Penelitian No
Penggunaan
Biaya (Rupiah)
1.
Penelusuran literatur
200.000,00
2.
Alat dan bahan penelitian
300.000,00
3.
Kolokium dan seminar hasil
500.000,00
4.
Penyusunan laporan
500.000,00 Total
Lampiran 4. Source Code dengan software Matlab Satu sel saraf function dsdt =subinagumo(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I = -0.35; x = s(1); y = s(2); % persamaan diferensialnya dsdt(1) = c*(x+y-x^3/3+I); dsdt(2) = -(x-a+b*y)/c;
Dua sel saraf terkopel function dsdt =subinagumocouple2(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3;
1.500.000,00
19
I =-0.45;%-1.45 Ts= -0.25; Vs= 2; k = 1; gs= 0.17; x1 y1 x2 y2
= = = =
s(1); s(2); s(3); s(4);
% persamaan diferensialnya dsdt(1) = c*(x1+y1-x1^3/3+I)-(x1-Vs)*gs*(1+exp(-k*(x2-Ts)))^(-1); dsdt(2) = -(x1-a+b*y1)/c; dsdt(3) = c*(x2+y2-x2^3/3+I)-(x2-Vs)*gs*(1+exp(-k*(x1-Ts)))^(-1); dsdt(4) = -(x2-a+b*y2)/c;
Tiga sel saraf terkopel function dsdt =subinagumocouple3(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I =-0.5;%-1.5 Ts= -0.25; Vs= 2; k = 1; gs=0.17;
x1 y1 x2 y2 x3 y3
= = = = = =
s(1); s(2); s(3); s(4); s(5); s(6);
% persamaan diferensialnya dsdt(1) = c*(x1+y1-x1^3/3+I)-(x1-Vs)*gs*((1+exp(-k*(x2-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)); dsdt(2) = -(x1-a+b*y1)/c; dsdt(3) = c*(x2+y2-x2^3/3+I)-(x2-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)); dsdt(4) = -(x2-a+b*y2)/c; dsdt(5) = c*(x3+y3-x3^3/3+I)-(x3-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)); dsdt(6) = -(x3-a+b*y3)/c;
Empat sel saraf terkopel function dsdt =subinagumocouple4(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I =-0.55;%-1.6 Ts= -0.25; Vs= 2; k = 1; gs=0.17; x1 y1 x2 y2 x3 y3
= = = = = =
s(1); s(2); s(3); s(4); s(5); s(6);
20
x4 = s(7); y4 = s(8); % persamaan diferensialnya dsdt(1) = c*(x1+y1-x1^3/3+I)-(x1-Vs)*gs*((1+exp(-k*(x2-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)); dsdt(2) = -(x1-a+b*y1)/c; dsdt(3) = c*(x2+y2-x2^3/3+I)-(x2-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)); dsdt(4) = -(x2-a+b*y2)/c; dsdt(5) = c*(x3+y3-x3^3/3+I)-(x3-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)); dsdt(6) = -(x3-a+b*y3)/c; dsdt(7) = c*(x4+y4-x4^3/3+I)-(x4-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x3-Ts)))^(-1)); dsdt(8) = -(x4-a+b*y4)/c;
Delapan sel saraf terkopel function dsdt =subinagumocouple8(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I =-0.81; Ts= -0.25; Vs= 2; k = 1; gs= 0.17; x1 = s(1); y1 = s(2); x2 = s(3); y2 = s(4); x3 = s(5); y3 = s(6); x4 = s(7); y4 = s(8); x5 = s(9); y5 = s(10); x6 = s(11); y6 = s(12); x7 = s(13); y7 = s(14); x8 = s(15); y8 = s(16); x9 = s(17); y9 = s(18); % persamaan diferensialnya dsdt(1) = c*(x1+y1-x1^3/3+I)-(x1-Vs)*gs*((1+exp(-k*(x2-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(2) = -(x1-a+b*y1)/c; dsdt(3) = c*(x2+y2-x2^3/3+I)-(x2-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(4) = -(x2-a+b*y2)/c; dsdt(5) = c*(x3+y3-x3^3/3+I)-(x3-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(6) = -(x3-a+b*y3)/c; dsdt(7) = c*(x4+y4-x4^3/3+I)-(x4-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x3-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(8) = -(x4-a+b*y4)/c; dsdt(9) = c*(x5+y5-x5^3/3+I)-(x5-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(10) = -(x5-a+b*y5)/c;
21
dsdt(11) = c*(x6+y6-x6^3/3+I)-(x6-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(12) = -(x6-a+b*y6)/c; dsdt(13) = c*(x7+y7-x7^3/3+I)-(x7-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)); dsdt(14) = -(x7-a+b*y7)/c; dsdt(15) = c*(x8+y8-x8^3/3+I)-(x8-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)); dsdt(16) = -(x8-a+b*y8)/c;
Enam belas sel saraf terkopel function dsdt =subinagumocouple16(t,s) dsdt = zeros(size(s)); % parameter a = 0.7; b = 0.8; c = 3; I =-1.4; Ts= -0.25; Vs= 2; k = 1; gs= 0.17; x1 = s(1); y1 = s(2); x2 = s(3); y2 = s(4); x3 = s(5); y3 = s(6); x4 = s(7); y4 = s(8); x5 = s(9); y5 = s(10); x6 = s(11); y6 = s(12); x7 = s(13); y7 = s(14); x8 = s(15); y8 = s(16); x9 = s(17); y9 = s(18); x10 = s(19); y10 = s(20); x11 = s(21); y11 = s(22); x12 = s(23); y12 = s(24); x13 = s(25); y13 = s(26); x14 = s(27); y14 = s(28); x15 = s(29); y15 = s(30); x16 = s(31); y16 = s(32); % persamaan diferensialnya dsdt(1) = c*(x1+y1-x1^3/3+I)-(x1-Vs)*gs*((1+exp(-k*(x2-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(2) = -(x1-a+b*y1)/c; dsdt(3) = c*(x2+y2-x2^3/3+I)-(x2-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x3Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(4) = -(x2-a+b*y2)/c;
22
dsdt(5) = c*(x3+y3-x3^3/3+I)-(x3-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(6) = -(x3-a+b*y3)/c; dsdt(7) = c*(x4+y4-x4^3/3+I)-(x4-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x3-Ts)))^(-1)+(1+exp(-k*(x5-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(8) = -(x4-a+b*y4)/c; dsdt(9) = c*(x5+y5-x5^3/3+I)-(x5-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x6Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(10) = -(x5-a+b*y5)/c; dsdt(11) = c*(x6+y6-x6^3/3+I)-(x6-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(12) = -(x6-a+b*y6)/c; dsdt(13) = c*(x7+y7-x7^3/3+I)-(x7-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x8-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(14) = -(x7-a+b*y7)/c; dsdt(15) = c*(x8+y8-x8^3/3+I)-(x8-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x9Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(16) = -(x8-a+b*y8)/c; dsdt(17) = c*(x9+y9-x9^3/3+I)-(x9-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(-k*(x2Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(18) = -(x9-a+b*y9)/c; dsdt(19) = c*(x10+y10-x10^3/3+I)-(x10-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x11-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(20) = -(x10-a+b*y10)/c; dsdt(21) = c*(x11+y11-x11^3/3+I)-(x11-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x12Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(22) = -(x11-a+b*y11)/c; dsdt(23) = c*(x12+y12-x12^3/3+I)-(x12-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(24) = -(x12-a+b*y12)/c; dsdt(25) = c*(x13+y13-x13^3/3+I)-(x13-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11-
23
Ts)))^(-1)+(1+exp(-k*(x12-Ts)))^(-1)+(1+exp(-k*(x14-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(26) = -(x13-a+b*y13)/c; dsdt(27) = c*(x14+y14-x14^3/3+I)-(x14-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11Ts)))^(-1)+(1+exp(-k*(x12-Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x15Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(28) = -(x14-a+b*y14)/c; dsdt(29) = c*(x15+y15-x15^3/3+I)-(x15-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11Ts)))^(-1)+(1+exp(-k*(x12-Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14Ts)))^(-1)+(1+exp(-k*(x16-Ts)))^(-1)); dsdt(30) = -(x15-a+b*y15)/c; dsdt(31) = c*(x16+y16-x16^3/3+I)-(x16-Vs)*gs*((1+exp(-k*(x1-Ts)))^(-1)+(1+exp(k*(x2-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x4-Ts)))^(-1)+(1+exp(-k*(x5Ts)))^(-1)+(1+exp(-k*(x6-Ts)))^(-1)+(1+exp(-k*(x7-Ts)))^(-1)+(1+exp(-k*(x8Ts)))^(-1)+(1+exp(-k*(x9-Ts)))^(-1)+(1+exp(-k*(x10-Ts)))^(-1)+(1+exp(-k*(x11Ts)))^(-1)+(1+exp(-k*(x12-Ts)))^(-1)+(1+exp(-k*(x13-Ts)))^(-1)+(1+exp(-k*(x14Ts)))^(-1)+(1+exp(-k*(x15-Ts)))^(-1)); dsdt(32) = -(x16-a+b*y16)/c;
