SIMULASI DINAMIKA SEL SARAF MENGGUNAKAN MODEL HINDMARSH-ROSE HINDMARSH ROSE
PUTRA WIRA KURNIAWAN
DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
INSTITUT PERTANIAN BOGOR BOGOR 2011
Putra Wira Kurniawan : Simulasi Dinamika Sel Saraf Menggunakan Model Hindmarsh-Rose. Dibimbing Oleh : Dr. Agus Kartono dan Dr. Ir. Irzaman M.Si.
ABSTRAK Dinamika penjalaran impuls pada sel saraf dalam tugas akhir ini adalah dinamika yang menggambarkan fenomena fisis yang dibangun dari persamaan Hindmarsh-Rose. Dinamika tersebut bersifat periodik sehingga memiliki titik-titik keseimbangan sistem dan bifurkasi lokal. Dengan perubahan input berupa arus searah maka didapatkan sebuah solusi numerik set persamaan differensial orde satu yang kemudian dapat dipecahkan secara analitik berupa matriks Jacobi untuk mendapatkan nilai eigennya. Metode RKF45 diimplementasikan pada perangkat lunak MATLAB yang berbasiskan GUI untuk menemukan solusi numerik dari persamaan differensial berorde yang terdapat dalam persamaan Hindmarsh-Rose. Keunggulan tersebut dapat lebih memudahkan pengguna secara umum untuk menganalisa dinamika sistem yang terjadi dalam sel saraf yang dibangun dari persamaan Hindmarh-Rose.
Kata kunci : Pengaruh berubahan arus searah I , penjalaran Impuls pada sel saraf, titik keseimbangan dan bifurkasi lokal, Graphical User Interface (GUI).
Judul Skripsi
: Simulasi Dinamika Sel Saraf dengan Menggunakan Model Hindmarsh-Rose
Nama
: Putra Wira Kurniawan
NIM
: G74070022
Menyetujui :
Pembimbing I,
Pembimbing II,
(Dr. Agus Kartono)
(Dr. Ir. Irzaman M.Si)
NIP : 197004211999031002
NIP : 196307081995121001
Mengetahui : Ketua Departemen,
(Dr. Ir. Irzaman M.Si) NIP : 196307081995121001
Tanggal Lulus :
SIMULASI DINAMIKA SEL SARAF DENGAN MENGGUNAKAN MODEL HINDMARSH-ROSE HINDMARSH ROSE
Skripsi Sebagai salah satu syarat untuk memperoleh gelar sarjana sains pada Departemen Fisika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Pertanian Bogor
Oleh:
PUTRA WIRA KURNIAWAN G74070022
DEPARTEMEN FISIKA FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
INSTITUT PERTANIAN BOGOR BOGOR 2011
DAFTAR RIWAYAT HIDUP PENULIS
Penulis dilahirkan di Jakarta Selatan pada tanggal 9 Februari 1989 oleh pasangan Ayah tercinta Mufrodi dan Ibunda tercinta Sri Yatmi. Penulis merupakan anak pertama dari dua bersaudara. Penulis menamatkan tingkat Taman Kanak-Kanak (TK) Seruni Pekanbaru pada tahun 1999 lalu melanjutkan di Sekolah Dasar (SD) 001 Rintis Pekanbaru dan menamatkannya pada tahun 2001. Penulis kemudian melanjutkan ke MTs Al-zaytun Indramayu Jawa Barat. Selama di tingkat tersebut penulis aktif mengikuti berbagai kegiatan ekstrakurikuler sebagai aktifis organisasi, olahragawan pencak silat, dan seniman nasyid. Setelah menamatkan jenjang MTs pada tahun 2004 penulis melanjutkan ke MA Al-zaytun Indramayu Jawa Barat. Pada tahun 2005 penulis pindah sekolah untuk menyelesaikan tingkat sekolah menengah atas di SMAN 1 Pekanbaru dan menamatkannya pada tahun 2007. Selama di tingkat menengah atas penulis juga aktif mengikuti kegiatan kerohanian dan mendapatkan perestasi belajar yang cukup memuaskan dengan masuk rangking 10 besar di kelas. Setelah tamat SMA, penulis melanjutkan studinya ke Institut Pertanian Bogor melalui jalur USMI dengan mengambil program studi Fisika. Selama kuliah penulis aktif menjadi anggota DPM TPB IPB (2007-2008), sekretaris I HIMAFI (Himpunan Mahasiswa Fisika) (2009-2010) dan aktif mengambil peran dalam berbagai kegiatan Himpro seperti menjadi staf Logstran Pesta Sains 2007, Tim Khusus Pesta Sains 2008 dan 2009, dan kegiatan departemen lainnya. Dalam keseharian penulis juga aktif dalam mendukung proses perkuliahan dan menularkan ilmunya dengan menjadi Asisten Praktikum Fisika Dasar (2008-2009) dan Asisten Fisika Komputasi (2011). Pada semester akhir ini penulis aktif mengamalkan ilmunya dengan menjadi pengajar di berbagai bimbingan belajar seperti PRIMAGAMA dan BKB Expert Multitalenta Bogor. Semua yang dilakukan penulis semata-mata ingin membekali diri agar layak untuk menggenggam masa depan yang lebih cerah dan semata-mata juga mengharapkan ridha Allah selaku hamba yang mencoba istiqomah dalam menjalankan syariat agama yaitu menuntut ilmu dan mengamalkannya. Moto hidup penulis adalah “Berbuatlah sesuatu didasari niat yang tulus dan ikhlas, lurus dan jangan pernah menyerah dengan banyak hal yang menyesatkan. Kesuksesan hanyalah milik orang yang mencoba dan selalu disertai doa, tawadlu dan keistiqomahan.”
i
KATA PENGANTAR
Puji syukur penulis penjatkan kehadirat Allah SWT, yang atas berkat dan rahmatNya maka penulis dapat menyelesaikan makalah usulan tugas akhir yang berjudul “Simulasi Dinamika Sel Saraf Menggunakan Model Hindmarsh-Rose”. Penulisan makalah hasil penelitian merupakan salah satu persyaratan untuk melakukan seminar hasil dan menyelesaikan penelitian tugas akhir. Penulis mohon maaf apabila dalam penulisan makalah ini masih banyak kekurangan-kekurangan baik teknis materi, mengingat akan kemampuan yang dimiliki penulis. Kritik dan saran dari semua pihak sangat penulis harapkan demi penyempurnaan dalam pembuatan makalah ini. Semoga makalah ini bermanfaat bagi semua. Kemudian penulis juga ingin mengucapkan terima kasih kepada orang tua yang selalu memberikan doa dan semangat. Kepada Dr. Agus Kartono dan Dr. Ir. Irzaman M.Si sebagai pembimbing skripsi dan pihak-pihak kelompok fisika bagian teori yang selalu memberikan bimbingan, saran, dan semangat untuk menyelesaikan penelitian ini. Kepada teman-teman fisika angkatan 44, 45, dan 43 beserta civitas akademika fisika lainnya yang saya banggakan. Akhir kata penulis berharap semoga Allah memberikan imbalan yang setimpal pada mereka yang telah memberikan bantuan, dan dapat menjadikan semua bantuan ini sebagai ibadah, Amiin Yaa Robbal ‘Alamiin.
Bogor, 07 Juni 2011
Penulis
ii
DAFTAR ISI DAFTAR TABEL....................................................................................................... iii DAFTAR GAMBAR .................................................................................................. iv DAFTAR LAMPIRAN............................................................................................... v BAB 1. PANDAHULUAN......................................................................................... 1 1.1 Latar Belakang............................................................................................ 1 1.2 Tujuan Penelitian........................................................................................ 1 1.3 Rumusan Masalah....................................................................................... 2 1.4 Hipotesis ..................................................................................................... 2 1.5 Keluaran...................................................................................................... 2 BAB 2. TINJAUAN PUSTAKA ................................................................................ 2 2.1 Morfologi Sel Saraf .................................................................................... 2 2.2 Fisiologi Sel Saraf....................................................................................... 2 2.3 Penjalaran Impuls Saraf.............................................................................. 5 2.4 Model Matematika Hinmarsh-Rose............................................................ 6 2.5 Dinamika Sistem dan Bifurkasi Lokal........................................................ 7 2.6 Metode Runge-Kutta-Fehlberg 45.............................................................. 9 BAB 3. BAHAN DAN METODE.............................................................................. 10 3.1 Tempat dan Waktu Penelitian..................................................................... 10 3.2 Peralatan ..................................................................................................... 10 3.3 Metode Penelitian ....................................................................................... 10 3.3.1 Analisis dinamika sistem....................................................................... 10 3.3.2 Pembuatan GUI..................................................................................... 10 BAB 4. HASIL DAN PEMBAHASAN ..................................................................... 11 4.1 Solusi Numerik Model H-R dengan RKF45............................................... 11 4.2 Analisis Bifurkasi dan Dinamika Sistem .................................................... 18 4.2 GUI untuk Dinamika Sistem H-R dengan MATLAB ................................ 22 BAB 5. KESIMPULAN DAN SARAN ..................................................................... 24 5.1 Kesimpulan ................................................................................................. 24 5.2 Saran ........................................................................................................... 25 DAFTAR PUSTAKA ................................................................................................. 26 LAMPIRAN................................................................................................................ 27
iii
DAFTAR TABEL Tabel 1.
Konsentrasi Ion di Dalam dan di Luar Sel Saraf Istirahat................. 3
Tabel 2.
Jenis-Jenis Titik Keseimbangan ........................................................ 8
Tabel 3.
Dinamika sistem penjalaran impuls pada sel saraf dengan variasi I input....................................................................... 16
Tabel 4.
Bifurkasi yang terjadi pada dinamika sistem dengan variasi I ................................................................................ 21
iv
DAFTAR GAMBAR Gambar 1.
Struktur sel saraf tunggal................................................................... 2
Gambar 2.
Distribusi tegangan listrik sepanjang membran sel saraf dalam kondisi istirahatnya................................................................. 3
Gambar 3.
Potensial aksi..................................................................................... 4
Gambar 4.
Penjalaran impuls saraf ..................................................................... 5
Gambar 5.
Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78A , a 3 , b 4, d 5 , 0.001 ,
1 2
dan xe Gambar 6.
5 1 ....................................................................... 15
Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78A , a 5 , b 14, d 5 , 0.001 , dan xe
Gambar 7.
1 2
5 1 ...................................................................... 15
Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78 , a 2 , b 9, d 5 , 0.001 , dan
xe
1 2
5 1
............................................................................... 16
Gambar 8.
Dinamika sistem penjalaran impuls pada sel saraf............................ 17
Gambar 9.
Analisa bifurkasi pada persamaan H-R ............................................. 20
Gambar 10.
Analisa bifurkasi dengan variasi I input ......................................... 22
Gambar 11.
GUI simulasi H-R model................................................................... 23
Gambar 12.
Penggunaan tombol ‘PLOT’ ............................................................. 23
Gambar 13.
Sifat dari konstanta xe sebagai titik keseimbangan bagi x0 pada saat I 0 pada system .............................................. 24
v
DAFTAR LAMPIRAN Lampiran 1.
Diagram Alir Penelitian .................................................................... 28
Lampiran 2.
Persamaan Hodgkin-Huxley untuk Potensial Aksi pada Satu Sel Saraf ........................................................................... 29
Lampiran 3.
Program Matlab untuk Metode RKF 45 dari persamaan Hindmarsh-Rose................................................................................ 30
Lampiran 4.
Program Matlab untuk GUI Fisik...................................................... 32
Lampiran 5.
Program Matlab untuk GUI Isi Program........................................... 36
1
BAB I PENDAHULUAN 1.1 Latar Belakang Sistem saraf merupakan salah satu sistem yang penting dalam mengatur fungsi kerja biologis. Sistem saraf didefinisikan sebagai suatu sistem koordinasi yang bertugas menyampaikan rangsangan dari reseptor untuk dideteksi dan direspon oleh tubuh.1 Sistem saraf memungkinkan makhluk hidup tanggap terhadap perubahan-perubahan yang terjadi di lingkungan luar maupun dalam. Pada hewan banyak ditemukan klasifikasi sistem saraf. Umumnya sistem saraf yang cukup terkenal adalah sistem saraf pusat dan tepi.1 Setelah menerima rangsangan atau stimulus baik yang berasal dari dalam maupun luar tubuh, rangsangan tersebut diteruskan ke sistem saraf pusat atau saraf tepi kemudian diintegrasikan dalam bentuk informasi guna menentukan respon yang akan diberikan oleh tubuh.1 Didalam sistem saraf terdapat sel saraf (neuron) yang merupakan bagian terkecil dalam suatu skema saraf dan berfungsi untuk menghantarkan informasi. Hampir seluruh jaringan mahluk hidup disusun oleh sel-sel saraf sebagai fungsi koordinasi dan pembawa sinyal neurotransmiter. Dalam memahami proses penjalaran impuls pada suatu sel saraf dibutuhkan pengetahuan mendasar tentang sifat konduktivitas membran dan mekanisme transportasi dalam sel saraf.1 Pada tugas akhir ini akan diteliti mengenai dinamika sistem dan simulasi penjalaran impuls pada sel saraf dengan menggunakan model Hindmarsh-Rose. Model Hindmarsh-Rose mengacu pada model analitik dan eksperimen yang telah dilakukan oleh A. L. Hodgkin dan A. F. Huxley (Hodgkin-Huxley), kemudian Hindmarsh-Rose mempelajari bagaimana penyederhanaan model matematika yang dilakukan oleh R.
Fitzhug dan J. Nagumo (FitzHugNagumo). Penyederhanaan tersebut mengambil persamaan gelombang Van der Pol dengan melakukan transformasi Leinard menjadi persamaan differensial dua variabel yang autonomous sehingga potensial aksi yang terbentuk merupakan activation dari saluran sodium (m) yang memiliki durasi yang sama dalam skala waktu tertentu, kemudian disaat inactivation saluran sodium (h) activation saluran potasium (n) sebagai fungsi deporlarisasnya.2,3 Mempelajari kelemahan dari model tersebut dimana menghasilkan potensial aksi yang sama perperiodik (inter-spike interval) dan tidak dapat menjelaskan hubungan fisis frekuensi dan arus.2 Hindmarsh-Rose mengubah fungsi linier yang ada pada persamaan Fitzhug-Nagumo menjadi persamaan kuadratik2, lalu dua tahun kemudian Hindmarsh-Rose menyempurnakan model analitiknya melalui persamaan ketiga agar model yang dibuat lebih menyerupai dinamika sistem yang ada pada penjalaran impuls di dalam sel saraf yakni bursting-firing adapting.3
1.2 Tujuan Penelitian a. Membuat simulasi dinamika sel saraf menggunakan model HindmarshRose dengan memvariasikan input berupa arus I searah sebagai sumber rangsangan. b. Menganalisa persamaan HindmarshRose yang merupakan persamaan diferensial terkopel pembentuk dinamika pada penjalaran impuls didalam sel saraf.
2
1.3 Perumusan Masalah a. Bagaimanakah pengaruh perubahan input berupa arus I searah sebagai sumber rangsangan terhadap penjalaran impuls didalam sel saraf dari persamaan Hindmarsh-Rose ? b. Bagaimanakah persamaan Hindmarsh-Rose membentuk dinamika dalam penjalaran impuls di dalam sel saraf ?
1.4 Hipotesis a. Input berupa arus I searah sebagai sumber rangsangan menunjukkan sejumlah fenomena penjalaran impuls spiking dan bursting pada sel saraf yang dibentuk dari persamaan Hindmarsh-Rose. b. Persamaan Hindmarsh-Rose merupakan persamaan differensial biasa yang bersifat autonomous dan saling terkopel yang dapat membentuk dinamika sistem didalam penjalaran impuls pada sel saraf.
terbungkus oleh suatu membran yang berisi cairan dengan nama aksoplasma.1 Pada perpanjangan ini terdapat selubung myelin, sel sachwann, dan Node of Renvier. Pada penjalaran impuls melalui sel saraf dendrit berfungsi menerima sinyal berupa rangsangan dan berfungsi sebagai sensor penerima dari sel saraf yang lainnya sedangkan akson berfungsi sebagai penghantar sinyal ke bagian sel saraf lain. Mekanisme penjalaran impuls melalui sel saraf terjadi sepanjang akson jika dan hanya jika rangsangan yang diterima oleh dendrit atau tubuh sel pada setiap waktu, intensitasnya berada pada ambang batas (threshold) atau lebih.1,5,7 Impuls akan mengalir dari tubuh sel menuju terminal akson. Sesampainya impuls saraf pada terminal akson, suatu subtansi saraf penghantar dilepaskan dan akan menyampaikan impuls ke penerima di sel berikutnya. Sambungan antara akson terminal dan dendrit sel saraf lain memungkinkan adanya sinapsis.5
1.5 Keluaran Keluaran yang diharapkan dari penelitian ini adalah sebuah simulasi berbasis Graphical User Interface (GUI) dari persamaan Hindmarsh-Rose yang dapat digunakan untuk menjelaskan dinamika sistem di dalam sel saraf.
BAB II TINJAUAN PUSTAKA 2.1 Morfologi Sel Saraf Sel saraf terdiri dari inti sel atau nukleus, tubuh sel atau sel somatik, diperpanjang oleh akson, terminal akson (ujung akson), dan sejumlah dendrit (lihat Gambar 1).1,5 Perpanjangan sel saraf dibentuk oleh akson yang merupakan salinan panjang tipis yang
©2003 Tippler, Inc.
Gambar 1. Struktur sel Saraf tunggal
2.2 Fisiologi Sel Saraf Seperti pada semua sel hidup, sel saraf memiliki kecenderungan mempertahankan kondisi intraseluler yang berbeda dengan lingkungan ekstraselularnya. Setiap sel saraf akan menghasilkan sejumlah ion negatif yang berada disekitar membran dalam sel dan ion positif yang berada mengintari bagian luar membran sel (Gambar 2a). Perbedaan besar muatan ion inilah menjadi dasar dalam menjelaskan mekanisme penjalaran impuls pada sel saraf. Selain itu sel syaraf juga memiliki sifat exitability (kemampuan merespon stimulus) dan conductivity (kemampuan menghantarkan sinyal).6
3
Tabel 1. Konsentrasi Ion di dalam dan di Luar Sel Saraf Istirahat.1 Jenis Ion +
©2003 Tippler, Inc.
(a) ©2007 Izhikevich, Inc.
(b) Gambar 2. (a) Potensial membrane dalam keadaan istirahat. (b) Distribusi tegangan listrik sepanjang membran sel saraf dalam kondisi istirahatnya.1,7 Sel saraf dalam kondisi tidak menghantarkan impuls berada pada keadaan istirahat atau (at rest). Keadaan ini dicirikan dengan adanya gradien konsentrasi ion-ion di bagian dalam dan di bagian luar membran pada nilai tertentu. Konsentrasi ion kalium (K+) dibagian dalam membran 30 kali lebih banyak daripada yang dibagian luar sedangkan konsentrasi ion Natrium (Na+) sepuluh kali lebih banyak dibagian luar membran sel dibandingkan dengan bagian dalam (Tabel 1). Adapun konsetrasi ion negatif (seperti Cl-) dan ion lainya tidak terdistribusi sempurna.1,7,8 Dalam mempertahankan kondisi ini sel saraf menggunakan mekanisme difusi pasif dan transportasi aktif. Ketidaksesuaian distribusi Na+ dan K+ terbentuk dari kebutuhan energi pemompaan Na+ K+, yang + + memindahkan Na dan K dari dalam dan luar sel (Gambar 2b).1,7
Na K+ ClLainnya
Konsentrasi (mmol/L) di dalam di luar 15 145 150 5 9 120 156 30
Saat sel saraf menghantarkan impuls, sejumlah protein pada membran sel akan berfungsi sebagai channel (saluran) yang memudahkan distribusi Na+ dan K+. Saluran Sodium (saluran Na+) dan Potassium (saluran K+) sangat spesifik untuk melewatkan ion-ion tertentu pada transportasi intra-membran selama proses penghantaran impuls sel saraf. Saat sel saraf berada pada keadaan istirahat, saluran Na+ bergantung pada tegangan tertutup, sehingga menjaga ketidaksamaan distribusi Na+. Tegangan tertutup dari saluran Na+ didapatkan dari gradien konsentrasi Na+ didalam dan diluar sel saraf dengan menggunakan persamaan Nernst (persamaan 1) maka didapatkan potensial membran (potensial tertutup) untuk Na+ 55mV (bergantung suhu lingkungan sel saraf untuk semua ion didalam sel saraf).9 Membran sel dalam keadaan istirahat tidak permeabel terhadap anion yang besar (atau terhadap jenis muatan negatif besar lainnya, seperti protein) dengan demikian kelebihan muatan negatif terbentuk tepat di bagian dalam permukaan membran sel. Dengan adanya konduktivitas dan polaritas dari membran sel maka terbentuklah suatu beda potensial antara bagian dalam dan luar membran sel yang disebut sebagai potensial istirahat (resting potensial). Nilai dari potensial istirahat tersebut sebesar -70mV (Gambar 2a).1,5,7,9
4
Potensial membran dalam keadaan istirahat tersebut dapat dihitung menggunakan persamaan Nerst sebagai berikut :
E
RT C i ln ZF C o
(1)
dengan :
E adalah potensial Nernst untuk ion x (diukur sebagai potensial membran atau potensial tertutup).
C o adalah konsentrasi dari ion C di luar sel.
Ci adalah konsentrasi dari ion C di dalam sel.
Z
adalah ion valensi dari C, untuk Z(K+) = 1 dan Z(Na+) = -1.
R adalah konstanta gas (8,314 J/molK).
T
istirahat, -70mV menjadi +30mV. Karena adanya gradien konsentrasi dan gradien listrik (Tabel 1 dan Gambar 2b), Na+ akan mengalir memasuki sel dan menimbulkan aliran arus listrik (persamaan 2).
I
q t
(2)
dengan I adalah arus yang mengalir Q adalah muatan ion t adalah waktu Fluks Na+ pada bagian dalam membran menghasilkan perubahan polaritas membran dan menyebabkan perubahan potensial listrik didalam membran menjadi lebih positif hingga mencapai +30 mV (menyimpang pada kisaran 100 mV dari potensial istirahat). Efek dari rangsangan ini diperlihatkan sebagai s2 pada Gambar 3a. ©2003 Tippler, Inc.
adalah suhu absolut.
F adalah konstanta Faraday (96485 C/mol). Potensial sel dalam keadaan istirahat dapat diganggu oleh rangsangan kimia dan fisis. Gangguan ini berpengaruh dalam waktu yang cukup singkat terhadap perubahan potensial istirahat, untuk setiap perubahan tersebut potensial membran sel membentuk suatu pola yang akan selalu kembali kepada potensial istirahat -70mV.1,10 Seperti yang dijelaskan sebelumnya pada semua rangsangan yang merambat dibutuhkan jika dan hanya jika rangsangan yang diterima oleh dendrit atau tubuh sel pada setiap waktu, intensitasnya berada pada ambang batas atau lebih. Rangsangan yang berada dibawah ambang batas tidak akan menyebabkan terbukanya saluran Na+ secara sempurna sehingga tidak menimbulkan depolarisasi (depolarization) maksimal dari potensial istirahat, hal ini diperlihatkan sebagai s1 pada Gambar 3a. Jika rangsangan yang diterima cukup, maka saluran Na+ akan membuka dengan sempurna serta menyebabkan depolarisasi dari potensial
(a)
©2006 Nathalie Corson, Inc.
(b) Gambar 3. (a) Potensial aksi. s1 adalah rangsangan dibawah nilai ambang, s2 rangsangan diatas ambang. (b) Mekanisme bukaan gerbang membran sodium dan potassium pada pembentukan potenisal aksi.1,5
5
Depolarisasi potensial membran diakibatkan muatan positif dari Na+ yang masuk kedalam membran sampai keadaan dimana konsentrasi Na+ berada pada keadaan seimbang di titik rangsangan semula. Saluran Potasium terbuka dan K+ keluar dari membran sel sehingga mengembalikan polaritas potensial membran ke keadaan semula (positif diluar dan negatif didalam). Keadaan ini mengkibatkan tertutupnya saluran Na+ dalam waktu yang cukup singkat dan sel saraf berada pada kondisi tidak dapat dirangsang lagi. Periode ini disebut sebagai periode pemulihan (repolarization). Transisi keadaan pada potensial membran ini disebut sebagai potensial aksi (spike). Selama masa pemulihan potensial istirahat dapat berada pada kondisi dibawah -70 mV keadaan ini disebut sebagai over pemulihan (hyperpolarization), keadaan ini dapat terjadi apabila saluran K+ membuka terlalu lama(Gambar 3b).2,5 Untuk mengembalikan konsentrasi ionion pada keadaan istirahat (setelah diberi stimulus) membran dapat melakukan transportasi aktif dengan memanfaakan tegangan tertutup dari masing-masing ion.
2.3 Penjalaran impuls saraf Pada fisiologi sel saraf telah dijelaskan bagaimana suatu pulsa listrik tunggal berupa rangsangan fisis berpengaruh pada tingkahlaku sel saraf. Namun penjalaran impuls saraf melalui sel saraf perlu dikaji lebih rinci lagi guna memahami mekanisme rangsangan untuk bisa sampai pada sistem saraf yang lebih kompleks.
©2003 Tippler, Inc.
(a) ©2003 Tippler, Inc.
(b) Gambar 4. (a) Penjalaran Impuls saraf, menyebabkan pertukaran K+ dan Na+ yang mengakibatkan timbulnya potensial aksi (b) di sepanjang akson.1 Setelah impuls yang diterima oleh dendrit atau badan sel melebihi batas ambang maka saluran Na+ terbuka untuk melewatkan Na+ memasuki membran sel sehingga menyebabkan depolarisasi lokal pada titik mula rangsangan, dan karena adanya gradien konsentrasi Na+ menyebabkan gerakan difusi-pasif ion yang yang berada pada daerah rangsangan (Gambar 4a). Karena adanya periode pemulihan (saluran K+ terbuka dan K+ keluar dengan membawa muatan positif), maka pada periode ini sebagian membran mengalami depolarisasi serta repolarisasi dan merambat pada satu arah tertentu saja, menjauhi tubuh sel saraf.1,5 Potensial aksi yang dihasilkan dari proses tersebut hanya terbentuk pada Renvier Node, potensial aksi melompat dengan cepat sepanjang akson (Gambar 4b) oleh adanya difusi ion-ion melalui aksoplasma dan cairan ekstraselular. Hal ini disebabkan aktivitas listrik pada sel saraf yang dilapisi myelin hanya terbatas pada Renvier Node daerah tersebut terdapat gradien konsentrasi yang cukup besar dari saluran ion yang bergantung pada tegangan tertutup. Impuls akan terus
6
bergerak hingga mencapai terminal, dan menyebabkan lepasnya neurotransmiter dari membran sel saraf.1 Neurotransmiter yang dihasilkan menjembatani rentang antar sel saraf sehingga seluruh proses dapat berulang atau diteruskan. Selama proses penghantaran impuls saraf, aliran listrik mengalir ke dalam dan keluar melalui membran tegak lurus terhadap arah perambatan impuls. Sehingga dengan asumsi ini seberapa jauhpun perpanjangan akson, impuls tidak pernah memerlukan penguatan, impuls akan terus merambat dengan kekuatan yang sama dari rangsangan awal.1 Sebagai tambahan diketahui bahwa sebagian akson diselimuti oleh sejumlah lapisan myelin yang terbentuk ketika sel-sel sachwann membungkus akson . Ruang antara lapisan myelin selebar 1 μm disebut Node of Renvier dan terbentuk pada setiap interval 1 sampai 2 mm sepanjang akson.1 Perambatan impuls melalui akson yang diselimuti lapisan myelin sedikit berbeda dengan perambatan melalui akson tanpa myelin. Lapisan myelin merupakan suatu insulator yang baik, sehingga ion tidak dapat mengalir menembus lapisan tersebut.1
2.4 Model Matematika Hindmarsh-Rose Dasar permodelan matematika dari sel saraf mengacu pada eksperimen yang telah dilakukan oleh Hodgkin-Huxley pada tahun 1952. Model Hodgkin-Huxley (H-H) menyatakan bahwa terdapat K+, Na+, dan ion lainnya yang mengalir melintasi membran serta perubahan nilai konduktivitas listrik membran terhadap ion-ion tersebut terhadap waktu sebagai fungsi dari potensial membran. Fisiologi tersebut dibuat dalam persamaan matematika dengan memformulasikan sebuah sistem persamaan differensial dengan empat
variabel yang potensial aksi.2,5,9
merepresentasikan
Pada kenyataannya pesamaan tersebut terdiri dari empat persamaan differensial terkopel nonlinier dengan enam fungsi dan tujuh konstanta 2,3, hal tersebut merupakan persamaan yang cukup kompleks untuk menemukan solusi numeriknya dan sangat sulit digunakan untuk simulasi dalam skala yang lebih kecil maka persamaan perlu disederhanakan. Pada tahun 1960 FitzHugh-Nagumo (FH-N) menyerderhanakan model H-H. FH-N mengubah sistem persamaan differensial empat variabel yang ada pada H-H menjadi dua persamaan differensial melalui persamaan gelombang Van der Pol dan transformasi Leinard. Kedua persamaan FH-N tersebut adalah : .
x a( y f ( x) I (t ))
(3a)
.
y b( g ( x) y)
(3b)
FH-N melakukan obsevasi penyederhanaan dengan menganggap V(t) merupakan activation saluran Na+ m(t) yang terbentuk saat potensial aksi, sementara itu inactivation saluran Na+ h(t) dan activation saluran n(t) berubah secara bersamaan dengan skala waktu yang sangat singkat. Variabel x menggambarkan potensial membran dan y adalah internal atau variabel pemulihan. Fungsi f(x) dan g(x) merupakan persamaan kubik dan linear, a dan b merupakan konstanta waktu dan I(t) adalah input arus yang bergantung waktu. Kemudian pada tahun 1982 Hindmarsh-Rose (H-R) menyederhanakan model H-H dengan mengambil konsep penyederhanaan yang telah dikenalkan oleh FitzHughNagumo (FH-N) sehingga bisa menjadi (slow-fast system). H-R mempelajari kelemahan dari model FH-N yang tidak dapat secara rinci menjelaskan fenomena (rapid firing)2,3 untuk itu H-R
7
merubah fungsi linier g(x) menjadi fungsi kuadratik (lihat persamaan 1).
f ( x) ax3 bx 2
(4a)
g ( x) c dx 2
(4b)
Meskipun model tersebut telah dapat menggambarkan (neural firing) secara jelas namun masih belum dapat menampilkan fitur biologi dari neuron seperti bursting dan adaptation. Dua tahun kemudian, H-R menambah persamaan ketiga dalam modelnya, agar dinamika dari model yang mereka buat lebih mendekati kondisi nyata. Persamaan terakhir yang mereka tambahkan dapat mengontrol jeda waktu diantara dua potensial aksi. .
x y ax2 x3 z I
(5a)
.
y 1 dx2 y
(5b)
.
z (b( x xe ) z )
(5c)
Dengan a = 3, b = 4 menjelaskan tentang pengaruh dari potensial membran saat (slow dynamics), d = 5 adalah konstanta yang ditentukan dari percobaan, adalah arus I 1.78 masukan, μ = 0.001 adalah recovery variable, xe
1 1 5 adalah titik 2
keseimbangan dari dua dimensi sistem H-R, disesuaikan dengan potensial ambang untuk memicu (bursting dan adaptation).5
2.5 Dinamika Sistem dan Bifurkasi Lokal Pada banyak kasus, dinamika sistem digunakan untuk memodelkan suatu fenomena yang memiliki hubungan sebab-akibat yang cukup kompleks. Dengan kata lain dinamika sistem dipergunakan untuk meningkatkan kemampuan pemahaman terhadap tingkah laku sistem yang dimunculkan dalam struktur tertentu.
Permodelan tersebut haruslah memiliki sifat dinamis (berubah terhadap waktu) dan struktur fenomenanya harus mengandung paling sedikit satu struktur umpan - balik (feedback structure). Oleh karena itulah model-model dinamika sistem diklasifikasikan ke dalam model matematik kausal (theory-like).10 Bifurkasi didefinisikan sebagai perubahan trayektori yang terjadi disekitar titik kritis, bifurkasi dicirikan dengan adanya perubahan jumlah titik kritis serta jenisnya akibat perubahan parameter yang terkandung di dalam suatu sistem persamaan.11 Titik kritis merupakan titik keseimbangan yang dimiliki oleh suatu dinamika sistem, titik-titik ini digunakan untuk menjelaskan bagaimana fenomena struktur dari dinamika sistem tersebut. Titik keseimbangan dapat ditentukan dengan mengenali persamaan differensial pembentuk dinamika sistem yang bersifat autonomous12, kemudian persamaan tersebut tidak berubah terhadap waktu atau sama dengan nol. Hal tersebut dapat ditinjau dari PDB berikut : .
x1 f1 ( x1 , x2 , ..., x N ) .
x2 f 2 ( x1 , x2 , ..., x N )
(6)
.
x N f N ( x1 , x2 , ..., x N ) Maka akan terdapat sejumlah titik keseimbangan xn xn ,0 yang diakibatkan secara f N ( x1,0 , ..., x n , 0 , ..., x N , 0 ) 0 serempak. Berdasarkan kenyataan tersebut sebuah titik keseimbangan dalam ruang–fasa terkait dengan solusi stationer sistem. Untuk mengetahui sifat dari titik kritis dapat dilakukan linierisasi dengan melakukan ekspansi Taylor terhadap f n di sekitar x n x n , 0 hingga orde pertama saja (dari ekspansi Taylor) :
8
f x n ( x n xn , 0 ) n xn n 1 .
N
....
(7)
Jenis Titik Keseimbangan
x xn , 0
Setelah mengetahui sifat dari titik keseimbangan tersebut dalam dinamika sistem dibutuhkan analisis lebih lanjut untuk mengetahui harga eigen. Harga eigen (λ) dibutuhkan untuk menentukan jenis bifurkasi yang terjadi dalam suatu sistem dinamik. Titik keseimbangan memberikan input dalam persamaan matriks jacobian J dari sistem :
f1 x1 f 2 J x1 f N x 1
f1 x2 f 2 x2 f1 x2
f1 x N f 2 x N f N x N
(8)
J adalah N N matriks yang diasumsikan sebagai matriks nonsingular yakni det J 0
(9)
maka dapat ditentukan harga eigen 1 , 2 , ..., N dengan menyelesaikan 12 persamaan berikut :
det( J I ) 0
Tabel 2. Jenis-Jenis Titik Keseimbangan12
(10)
Analisa dinamika sistem melibatkan sejumlah peran titik keseimbangan trayektori dan perubahannya (bifurkasi). Kedua analisa tersebut memiliki sifat dan jenis tertentu untuk menggambarkan trayektori (lintasan) sistem sebelum dan sesudah bifurkasi. Karakteristik titik keseimbangan yang dapat diketahui setelah menemukan harga eigen seperti tertera pada tabel berikut :
Harga Eigen
Titik Node
rill, untuk kasus λ1, λ2 > 0 (atraktor negatif) menjauhi titik keseimbangan, untuk kasus λ1, λ2 < 0 (atraktor positif) menuju titik keseimbangan.
Titik Sadel
rill, λ1 > 0, λ2 < 0 atau sebaliknya.
Titik Center
imajiner.
Titik Fokus
bilangan kompleks (terdiri dari rill dan imajiner) untuk kasus μ > 0 atraktor negatif menjauhi titik kritis, untuk kasus μ < 0 atraktor positif menuju titik kritis.
Jenis dari bifurkasi lokal pada dinamika sistem cukup banyak dikenali diantaranya bifurkasi Sadel-Node, bifurkasi Trans-Kritikal, bifurkasi PitchFork, dan bifurkasi Poincare-AndronovHopf (bifurkasi Hopf).12,13,14 Namun tidak semua jenis bifurkasi tersebut terbentuk dan menjadi fenomena struktur dari suatu dinamika sistem. Dinamika sistem yang linier merupakan kondisi stabil yang berpusat pada titik keseimbangan, namun pada kenyataannya dinamika sistem di alam menunjukkan tingkah laku yang tidak harmonis atau non-linier. Tingkah laku tersebut dapat membentuk suatu struktur fenomena yang teratur ditengah ketidak teraturan atau sewaktu-waktu dapat terlihat tidak teratur ditengah keteraturan (teori chaotik).12,13 Teori chaotik meninjau sifat ketidakstabilan pola osilasi sistem yang bergantung pada prediksi waktu. Lebih umumnya osilasi ini dapat berbentuk suatu osilasi harmonis pada selang tertentu dan dapat berbentuk tidak harmonis pada selang waktu yang lainnya. Teori yang mendasari dari osilasi chaotik meninjau suatu keadaan harmonis menuju ke keadaaan tidak harmonis atau sebaliknya. Sistem dapat dikatakan chaotik bila sistem tersebut dapat ditentukan secara deteministik, memiliki tingkah laku tidak periodik dalam masa yang cukup lama, dan menampakkan
9
perubahan yang bergantung pada pengaturan kondisi awal dari sistem (dalam hal ini berupa input atau inisialisasi kondisi sistem).14
2.6 Metode Runge-Kutta-Fehlberg 45 (RKF45) Metode Runge-Kutta-Fehlberg 45 (RKF45) merupakan salah satu metode untuk memecahkan solusi numerik dari persamaan differensial biasa berorde tinggi. Metode ini membuat suatu pencacahan dengan menggunakan selang h dengan enam konstanta yang akan diaproksimasikan kedalam persamaan orde 4 dan 5 sehingga lebih teliti dibanding metode lainnya.17 Keenam konstanta berperan untuk memprediksi harga solusi yang diinginkan pada dua keadaan atau lebih dari suatu sistem sedemikian rupa sehingga galat pembulatan dapat diminimisasi sampai orde 4 hingga orde 5.15,16 Metode ini ditampilkan dalam persamaan berikut :
k1 hf t k , yk ,
(11a)
1 1 k 2 hf t k h, y k k1 , 4 4
(11b)
3 3 t k 8 h, y k 32 k1 , (11c) k 3 hf 9 k 32 2
12 t k 13 h, y k 1932 k 2197 1 , k 4 hf 7200 k 2 2197 7296 k 3 2197
(11d)
439 t k h, y k 216 k1 3680 k 5 hf 8k 2 k 3 , (11e) 513 845 4104 k 4 1 8 t k 2 h, y k 27 k1 3544 k 6 hf 2k 2 k3 , (11f) 2565 11 1859 4104 k 4 40 k 5 f merupakan fungsi yang dapat dibentuk dari suatu persamaan differensial. Konstanta-konstanta tersebut kemudian dapat disubstitusikan ke dalam persamaan orde 4 (persamaan 12) dan orde 5 (persamaan 13) sebagai berikut :
y k 1
12 1408 y k k k 216 1 2565 3 2197 k 1 k 4101 4 5 5
(12)
k 6 yang tidak substitusikan dalam orde 4 namun disubstitusikan pada orde 515
16 6656 y k 135 k1 12825 k 3 9 28516 y k 1 k 4 k5 (13) 56430 50 2 55 k 6 Untuk k 2 tidak disubtitusikan dalam kedua persamaan orde tersebut hanya digunakan sebagai determinasi untuk konstanta lainnya.15 Metode ini diaplikasikan pada perangkat lunak MATLAB untuk membantu menyelesaikan solusi numerik yang kompleks dan simulasi berbasis Graphical User Interface (GUI). GUI MATLAB dibangun dari fisik dasar yaitu ‘figure’ sebagai alokasi
10
media utama yang kemudian dilanjutkan dengan ‘uicontrol’ sebagai tampilan dalam menjalankan fungsi khusus dan spesifik atau sebagai fungsi unit interface-nya, misalnya fungsi pushbutton untuk membentuk tombol yang dapat di tekan dan fungsi edit untuk membentuk kolom isian. Fungsi untuk menampilkan keluaran dari GUI adalah fungsi axes.18 MATLAB dikenal sebagai perangkat lunak yang digunakan untuk komputasi teknik dan sains. Perangkat lunak ini bekerja pada operasi dasar matematika (kalkulator) dan algoritma pemograman tertentu. Sehingga MATLAB dapat beroperasi untuk menemukan solusi numerik baik berupa skalar, vektor, maupun pengolahan data lebih lanjut.18
BAB III BAHAN DAN METODE 3.1 Tempat dan Waktu Penelitian Penelitian ini dilakukan di Laboratorium Fisika Teori dan Komputasi, Departemen Fisika, Fakultas Matematika dan Ilmu Pengetahuan Alam, Institut Pertanian Bogor dari bulan September 2010 sampai dengan April 2011.
3.2 Peralatan Peralatan yang digunakan adalah sebuah laptop dengan processor Intel Core2Duo T6400 2.0 GHz / 800 Mhz FSB, L2 Chace / 45 nm technology, 2 GB RAM. Perangkat lunak yang digunakan adalah MS. Office 2007 dan MATLAB R2009b. Sebagai pendukung penulis menggunakan sumber literatur, yaitu jurnal-jurnal ilmiah, buku-buku, dan sumber-sumber lain yang terkait.
3.3 Metode Penelitian Penelitian ini dimulai dengan membuat simulasi penjalaran impuls pada sel saraf menggunakan model Hindmarsh-Rose (lihat 2.4 Model Matematika Hindmarsh-Rose). Simulasi model ini menggunakan metode numerik Runge-Kutta-Fehlberg 45 (lihat 2.6 Runge-Kutta-Fehlberg 45). Kemudian dilakukan analisa terhadap persamaan pembentuk dinamika sistem (lihat 2.5 Dinamika Sistem dan Bifurkasi Lokal) yang terjadi pada persamaan Hindmarsh-Rose.
3.3.1 Analisa dinamika sistem Persamaan Hindmarsh-Rose yang terdiri dari persamaan differensial biasa dengan tiga variabel yang autonomous atau eksplisit akan dianalisa dengan mencari titik keseimbangannya menggunakan (persamaan 5 hal.7). Setelah itu dengan menggunakan ekspansi Taylor (persamaan 7 hal.8) dan penemuan harga eigen (persamaan 8 hal.8) untuk mengetahui bifurkasi lokal yang terjadi dalam persamaan tersebut. Bifurkasi lokal yang terjadi dalam dinamika sistem penjalaran impuls pada sel saraf di ketahui dengan meninjau perubahan dan jenis dari titik keseimbangan. Variasi inputan berupa besarnya arus (impuls) akan mempengerahui dinamika sistem dan bifurkasi lokal dari persamaan tersebut dalam membentuk penjalaran impuls pada sel saraf.
3.3.2 Pembuatan GUI Pembuatan GUI dilakukan dengan software MATLAB R2009b dengan tujuan mempermudah menganalisa dinamika sistem yang terjadi pada persamaan Hindmasrh-Rose. GUI yang dibentuk memungkinkan untuk dilakukan perubahan konstanta dan input sehingga setiap perubahan tersebut dapat menampilkan struktur fenomena yang berbeda.
11
BAB IV HASIL DAN PEMBAHASAN 4.1 Solusi Numerik Model H-R dengan RKF45 Model H-R yang terbentuk dari tiga persamaan differensial orde satu yang saling berhubungan atau terkopel. Persamaan tersebut bersifat autonomous yang berarti berdiri sendiri secara eksplisit dari variabel penyusunnya tanpa variabel waktu yang mempengaruhinya, meskipun diturunkan dari fungsi waktu. H-R mengembangkan persamaan analitiknya dari tahun 1982 s.d. 1984 untuk meniru dinamika sistem yang terjadi dalam penjalaran impuls pada sel saraf. Model ini merupakan model penyederhanaan (simplify model) yang mengacu pada fenomena eksperimen yang telah dilakukan oleh Hodgkin-Huxley atau yang lebih dikenal dengan model H-H.2,3,4,5 Penjaralan impuls saraf pada satu sel saraf memperlihatkan dinamika yang kompleks dan khas. Fenomena tersebut diantaranya adalah adanya potensial istirahat, potensial ambang (threshold), hingga timbulnya potensial aksi atau spike.1 Namun hal tersebut hanya merupakan sebagian kecil dari fenomena-fenomena lain yang lebih kompleks dan menarik untuk dipelajari. Pada mulanya H-R menemukan persamaan penjalaran impuls saraf dengan dua buah variabel peubah dari persamaan differensial biasa orde satu yang autonomous, persamaan ini merupakan model H-R yang lebih dikenal dengan slow fast system.5 Sistem persamaan ini dianggap berhasil untuk meniru fenomena spike penjalaran impuls pada sel saraf dan menutupi kelemahan dari model penyederhanaan yang ditemukan oleh FitzHug-Nagumo (FH-N). Namun model ini masih memiliki kekurangan untuk menjelaskan
fenomena firing-adapting. Sehingga pada tahun 1984 H-R memformulasikan persamaannya menjadi tiga persamaan differensial orde satu yang saling terkopel. Penambahan satu peubah disamping dua peubah sebelumnya dimaksudkan untuk memperlihatkan tingkah laku penjalaran yang lebih nyata seperti (bursting dan adaptation).3 Persamaan analitik yang dibangun oleh H-R merupakan persamaan nonlinier yang tidak mudah dipecahkah untuk mencari solusi fisisnya, sehingga didalam ilmu sains penyelesaian solusi tersebut dapat dilakukan secara numerik.9 Banyak metode numerik yang sudah dikenal dapat memecahkan solusi dari persamaan differensial biasa, diantaranya metode Euler, Hund, dan Runge-Kutta. Dalam kasus ini penyelesaian numerik dari model H-R dilakukan dengan metode RKF45 (lihat 2.6 Runge-Kutta-Fehlberg 45). Metode ini adalah penyempurnaan dari metode Runge-Kutta orde 4 dan teruji lebih baik dengan menghasilkan nilai kesalahan (error value) yang paling kecil dari metode sebelumnya. Implementasinya adalah memasukkan persamaan H-R ke dalam metode RKF45. Metode RKF45 mencacah sejumlah titik sehingga menghasilkan trayektori dalam ruang fasa tertentu berdasarkan persamaan berupa fungsi yang dibentuk dari model H-R. Terdapat tiga fungsi dari persamaan H-R yang bisa dibentuk dari persamaannya (persamaan 5 hal. 7):
f x y ax 2 x 3 z I
(14a)
f y 1 dx 2 y
(14b)
f z (b( x xe ) z )
(14c)
Ketiga fungsi yang telah dibentuk tersebut disubstitusikan ke metode RKF45 sebagai berikut (lihat persamaan 11 dan 13 pada hal. 9):
12
kx1 hf x xk , y k , z k ,
(15a)
ky1 hf y x k , y k , z k ,
(15b)
kz1 hf z xk , yk , z k , 1 1 1 kx2 hf x x k kx1 , y k ky1 , z k kz1 , 4 4 4
(15c)
1 1 1 ky 2 hf y x k kx1 , y k ky1 , z k kz1 , 4 4 4
(15d)
1 1 1 kz 2 hf z x k kx1 , y k ky1 , z k kz1 , 4 4 4
(15e)
3 9 3 9 kx1 kx2 , y k ky1 ky 2 , xk 32 32 32 32 , kx3 hf x 3 9 kz1 kz 2 zk 32 32
(15f)
3 9 3 9 kx1 kx2 , y k ky1 ky 2 , xk 32 32 32 32 , ky3 hf y 3 9 kz1 kz 2 zk 32 32
(15g)
3 9 3 9 kx1 kx2 , y k ky1 ky 2 , xk 32 32 32 32 , kz3 hf z 3 9 kz1 kz 2 zk 32 32
(15h)
1932 7200 7296 kx1 kx2 kx3 , xk 2197 2197 2197 1932 7200 7296 kx4 hf x y k ky1 ky 2 ky3 , , 2197 2197 2197 z 1932 kz 7200 kz 7296 kz 1 2 3 k 2197 2197 2197
(15i)
1932 7200 7296 kx1 kx2 kx3 , xk 2197 2197 2197 1932 7200 7296 ky 4 hf y y k ky1 ky 2 ky3 , , 2197 2197 2197 z 1932 kz 7200 kz 7296 kz 1 2 3 k 2197 2197 2197
(15j)
13
1932 7200 7296 kx1 kx2 kx3 , xk 2197 2197 2197 1932 7200 7296 ky 4 hf Z y k ky1 ky 2 ky3 , , 2197 2197 2197 z 1932 kz 7200 kz 7296 kz 1 2 3 k 2197 2197 2197
(15k)
439 3680 845 kx1 8kx 2 kx3 kx 4 , xk 216 513 4104 439 3680 845 kx5 hf x y k ky1 8ky 2 ky 3 ky 4 , , 216 513 4104 439 3680 845 z kz1 8kz 2 kz 3 kz 4 k 216 513 4104
(15l)
439 3680 845 kx1 8kx2 kx3 kx 4 , xk 216 513 4104 439 3680 845 ky5 hf y y k ky1 8ky 2 ky3 ky 4 , , 216 513 4104 z 439 kz 8kz 3680 kz 845 kz 1 2 3 4 k 216 513 4104
(15m)
439 3680 845 kx1 8kx2 kx3 kx4 , xk 216 513 4104 439 3680 845 kz 5 hf z y k ky1 8ky 2 ky3 ky 4 , , 216 513 4104 z 439 kz 8kz 3680 kz 845 kz 1 2 3 4 k 216 513 4104
(15n)
8 3544 1859 11 kx1 2kx2 kx3 kx4 kx5 , xk 27 2565 4104 40 8 3544 1859 11 kx6 hf x yk ky1 2ky2 ky3 ky4 ky5 , , 27 2565 4104 40 8 3544 1859 11 z kz1 2kz 2 kz3 kz 4 kz5 k 27 2565 4104 40
(15o)
8 3544 1859 11 kx1 2kx2 kx3 kx4 kx5 , xk 27 2565 4104 40 8 3544 1859 11 ky6 hf y y k ky1 2ky2 ky3 ky4 ky5 , , 27 2565 4104 40 8 3544 1859 11 z kz1 2kz2 kz3 kz4 kz5 k 27 2565 4104 40
(15p)
14
8 3544 1859 11 kx1 2kx2 kx3 kx4 kx5 , xk 27 2565 4104 40 8 3544 1859 11 kz6 hf z y k ky1 2ky2 ky3 ky4 ky5 , , 27 2565 4104 40 8 3544 1859 11 z kz1 2kz2 kz3 kz4 kz5 k 27 2565 4104 40 Konstanta-konstanta tersebut dibangun berdasarkan metode RKF45 yang akan membentuk perubahan titik-titik untuk trayektori dengan menggunakan persamaan orde 5 (lihat 2.7 Metode Runge-Kutta-Fehlberg 45) disajikan sebagai berikut:
xk 1
16 xk 135 kx1 6656 kx 12825 3 28561 kx 56430 4 9 2 kx5 kx6 55 50
16 y k 135 ky1 6656 ky 12825 3 y k 1 28561 ky 56430 4 9 2 ky5 ky6 55 50
z k 1
16 z k 135 kz1 6656 kz 12825 3 28561 kz 56430 4 9 2 kz5 kz6 55 50
(16a)
(16b)
(16c)
Untuk k 1 dimana t1 , x1 , y1 , z1 nilai yang disesuaikan sebagai inisial dan n (jumlah iterasi) merupakan batas iterasi t k 1 t k h k ( k 1 s.d. n ).
(15q)
merupakan perubahan waktu dengan
h
t , n
t adalah
selang
waktu
diinginkan untuk meninjau trayektori dari persamaan H-R. Perlu diketahui dalam metode RKF45 syarat untuk meninjau iterasi titik data dengan akurasi tinggi diperlukan nilai n yang jauh lebih besar dari t . Penjelasan diatas merupakan tahapan untuk menemukan sebuah solusi numerik berupa grafik trayektori untuk membuat simulasi sistem dinamika penjalaran impuls di dalam sel saraf dengan menggunakan perangkat lunak MATLAB yang dapat dilihat pada Gambar 5. Grafik tersebut merupakan simulasi yang memperlihatkan fenomena terbentuknya potensial aksi dengan nilai 3 satuan potensial sehingga dapat disesuaikan dengan nilai yang sebenarnya dari nilai potensial aksi berdasarkan eksperimen H-H. Konstanta-konstanta yang terdapat dalam persamaan H-R dapat dijelaskan sebagai berikut (lihat 2.5 Model Matematika Hindmarsh-Rose) : a dan b yang menjelaskan tentang pengaruh dari potensial membran saat (slow dynamics). Slow dynamics adalah keadaan saat terbukannya saluran Na+ yakni saat potensial didalam membran sel saraf mulai meningkat menuju potensial aksi. Perubahan variabel ini akan merubah nilai potensial aksi dan lamanya proses perubahan potensial menuju potensial aksi. Dalam kasus ini perubahan nilai akan membuat perbandingan baru untuk menerka kuantitas dari proses pembentukan potensial aksi, bahkan dengan perubah sampai tertentu harga tertentu dapat membuat trayektori simulasi potensial
15
aksi yang tidak sesuai dengan fenomena yang ada. Aplikasinya ada pada keadaan sel saraf mahluk hidup yang memiliki fenomena berbeda untuk menyatakan kuantitas proses terbentuknya potensial aksi. Perubahan varibel a dan b dibandingkan antara Gambar 5, Gambar 6, dan Gambar 7. Perubahan variabel tersebut berkaitan dengan karakteristik dari bukaan pintu Na+ setiap sel saraf dan bergantung juga pada kondisi fisik dari sel saraf. Semakin lama bukaan dari pintu Na+ maka akan mengakibatkan pembentukan potensial aksi yang semakin terlambat pada Gambar 7. Variabel d adalah konstanta yang ditentukan dari percobaan sehingga tidak ada penjelasan khusus dari perubahan variabel ini pada terbentuknya potensial aksi, namun apabila variabel ini diubah maka trayektori yang terbentuk akan menunjukkan fenomena lain diluar fenomena yang dicari dari simulasi ini. I adalah arus masukan, dalam hal ini bukan merupakan besarnya arus yang persis diberikan pada suatu sel saraf hidup, namun merupakan arus dugaan dari kuantitas fenomena potensial aksi yang telah didapatkan sebelumnya. Variabel adalah recovery variable dan xe adalah titik keseimbangan dari dua dimensi sistem H-R yang akan bernilai dirinya sendiri apabila tidak ada inputan didalam sistem (dapat dilihat dengan menggunakan GUI pada sub-bab 4.3). 5
Gambar 5. Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78A , a 3 ,
d 5,
b 4,
xe
1 2
0.001 ,
dan
5 1 .
Potensial aksi merambat melalui membran sel saraf dan membentuk potensial aksi.
Gambar 6. Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78 A , a 5 ,
d 5,
b 14,
xe
1 2
0.001 ,
dan
5 1 .
Perubahan potensial aksi saat slow dynamics akibat perubahan nilai a dan b .
16
Gambar 7. Simulasi potensial aksi dengan perangkat lunak MATLAB dari persamaan H-R. I 1.78 A , a 2 ,
d 5,
b 9,
xe
1 2
0.001 ,
dan
5 1 .
Pada variasi nilai tertentu perubahan nilai a dan b membuat simulasi diluar fenomena pada I yang sama.
Dinamika spiking merupakan pola potensial aksi yang menjalar sepanjang akson sel saraf tanpa mengalami bursting yang artinya tidak ada pengulangan secara serempak pada suatu kumpulan perambatan potensial aksi. Perilaku bursting dianggap pola penjalaran yang khas, karena menumpuknya potensial aksi dalam satu waktu penjalaran yang singkat, dimisalkan seperti laju tranfer sinyal neurotransmiter persatuan panjang akson. Adapun perilaku fast spiking merupakan perilaku umum pada sel saraf yang menunjukkan pengulangan spiking secara normal dalam periode yang lebih cepat. Fast Spiking terjadi pada sel saraf dengan hantaran impuls searah. Dengan memvariasikan I maka didapatkan sejumlah data yang tersaji pada Tabel 3 berikut :
Tabel 3. Dinamika sistem penjalaran impuls pada sel saraf dengan variasi I input Selang variasi I
0 I 1.33 I 1.34 1.35 I 1.41 1.42 I 1.58 1.59 I 1.77 1.78 I 1.95 1.96 I 2.13 2.14 I 2.28 2.29 I 2.42 2.43 I 2.55 2.56 I 2.66 2.67 I 2.76 2.77 I 2.85 2.85 I 2.93 2.93 I 3.00 3.01 I 3.07 3.08 I 3.13 3.14 I 3.18 3.19 I 3.24 3.25 I 3.27 I 3.28 3.29 I 3.30 I 3.31 3.32 I
Dinamika sistem yang terjadi Spiking Bursting dengan 4 spike Bursting dengan 5 spike Bursting dengan 6 spike Bursting dengan 7 spike Bursting dengan 8 spike Bursting dengan 9 spike Bursting dengan 10 spike Bursting dengan 11 spike Bursting dengan 12 spike Bursting dengan 13 spike Bursting dengan 14 spike Bursting dengan 15 spike Bursting dengan 16 spike Bursting dengan 17 spike Bursting dengan 18 spike Bursting dengan 19 spike Bursting dengan 20 spike Bursting dengan 21 spike Bursting dengan 22 spike Bursting dengan 25 spike Bursting dengan 23 spike Chaotic Bursting Spiking
17
(a)
(c)
(b)
(d)
Gambar 8. Dinamika sistem penjalaran impuls pada sel saraf. (a) Untuk I 1.34A Busting dengan 4 spike. (b) Untuk I 3.08A Busting dengan 19 spike. (c) Untuk
I 3.31A Chaotic Busting. (d) Untuk I 3.60A Fast Spiking.
18
Tabel diatas merupakan sebuah hasil simulasi yang memperlihatkan fenomena-fenomena terbentuknya potensial aksi didalam penjalaran impuls pada sel saraf. Terlihat dari tabel semakin besarnya inputan arus searah I mengakibatkan semakin terbentuknya fenomena bursting yang dimulai dari fenomena dasar yakni spiking. Terlihat 0.05A pada kenaikan selang terbentuklah sebuah spike tambahan I 3.27A tepat hingga pada berakhirnya penjumlahan spike pada bursting secara linear. Perlu diketahui fenomena bursting merupakan sebuah fenomena terbentuknya penumpukan jumlah potensial aksi atau spike pada waktu tertentu dan membentuk pola yang teratur sepanjang penjalaran impuls didalam sel saraf. Fenomena ini sangat umum terjadi dan menjadi hal yang penting untuk memperlihatkan kondisi normal dan tidaknya suatu sel saraf dari karakteristiknya menyampaikan sinyal neutransmiter 5. Jumlah spike yang terbentuk bergantung dari besarnya arus yang di berikan pada setiap jenis dari sel saraf. Analisa penambahan linier arus dan pembentukan spike bursting terhenti saat arus berada I 3.27A karena pada
pada pada pada saat
I 3.28A potensial aksi menjalar dengan pola spike bursting yang meningkat menjadi 25 spike dan kembali menurun menjadi 23 spike fenomena selanjutnya merupakan fenomena chaotic bursting yang berarti terjadi fenomena kacau pada bursting. Terjadinya hal tersebut merupakan temuan yang terkait dengan keterbatasan persamaan H-R dalam menjelaskan fenomena penambahan jumlah spike yang linier untuk setiap penambahan jumlah arus. Namun hal tersebut tidak sepenuhnya menjadi kelemahan dari model ini, pada kenyataannya model ini juga menjelaskan tentang kekacauan pulsa dari potensial aksi yang menjalar dalam sebuah sel saraf akibat berlebihnya impuls yang diberikan 1,
sehingga mengakibatkan sel saraf mengalami chaos 9 Sementara itu untuk I 3.32A simulasi persamaan memperlihatkan fenomena spiking secara periodik yang akan semakin cepat dengan penambahan arus.
4.2 Analisa Bifurkasi dan Dinamika Sistem Tingkah laku potensial aksi (singgle spiking) yang telah dibahas sebelumnya merupakan sebagian kecil dari fenomena penjalaran impuls pada sel saraf. Masih banyak fenomena lainnya yang terjadi dan merupakan fenomena yang wajib untuk diketahui dan dipelajari, diantaranya adalah tingkah laku bursting, chaotic bursting, bursting adapting, dan lain-lain. Tingkah laku tersebut dapat dilihat dari grafik simulasi dengan variasi inputan I arus pada persamaan atau impuls (rangsangan) dan keadaan dari sel saraf (konstanta dalam persamaan) sendiri. Dalam kasus ini analisis bifurkasi dan dinamika sistem merupakan analisis yang yang digunakan untuk menganalisa persamaan sehingga dapat membentuk beberapa fenomena diatas. Seperti yang telah disebutkan sebelumnya pengertian bifurkasi adalah perubahan trayektori dalam ruang fasa (dalam hal ini potensial aksi) yang terjadi disekitar titik kestabilan. Titik kestabilan merupakan titik yang ditentukan dengan membuat semua fungsi dari persamaan differensial terkopel sama dengan nol. Kasus ini menggunakan tiga persamaan differensial H-R 11:
f x y ax 2 x 3 z I 2 f y 1 dx y f z bx xe z fx 0 f y 0 fx 0
(17a)
(17b)
19
y ax 2 x 3 z I 0 1 dx 2 y 0 b x x z 0 e
(17c)
Kemudian didapatkan persamaan polinomial sebagai berikut :
x 3 d a x 2 bx bx e I 1 0
(18)
untuk I 1.78 , a 3 , b 4, d 5 ,
1 0.001 , dan xe 2
5 1 ,
Akar-akar dari persamaan tersebut adalah titik kestabilan x0 dalam ruang fasa yang bergantung pada variasi I . Dengan mengetahui titik kestabilan untuk x maka didapatkan titik kestabilan untuk y dan z yakni :
x0 y0 z 0
- 1.2125 - 6.3504 1.6223
Tingkah laku penjalaran impuls pada sel saraf untuk setiap besar inputan (variasi I ) akan menunjukan nilai titik kestabilan yang selalu berubah, hal tersebut berakibat pada timbulnya polapola penajalaran impuls yang bervariasi dalam selang waktu tertentu dan memastikan adanya fenomena-fenomena yang telah disebutkan sebelumnya (bursting, chaotic bursting, bursting adapting, dan lain-lain) pada suatu harga I dalam selang tertentu pula. Sehingga dapat dianalisa dan diprediksi bagaimana dinamika sistem penjalaran impuls yang akan terjadi di dalam sel saraf. Namun sebelum menganalisa dinamika sistem lebih lanjut perlu dilakukan pengkarakterisasian dari harga-harga keseimbangan untuk mengetahui jenis bifurkasi dari pembentuk dinamika sistemnya. Harga eigen merupakan harga yang didapatkan dari titik kritis pada suatu sistem dengan melakukan substitusi nilai
salah satu nilai atau seluruhnya dari x0 , y0 , z0 kedalam suatu matriks untuk melakukan pengkarakterisasian. Karakterisasi disini adalah untuk mengetahui jenis bifurkasi apa yang telah dialami suatu sistem dari persamaan eksak pembentuknya. Harga eigen karakteristik didapatkan dengan melakukan linearisasi persamaan H-R untuk membentuk Jacobian Matrix (lihat persamaan 8 hal. 8) 11 :
3 x 2 2ax 1 1 1 0 , J 2dx b 0
(19a)
persamaan matriks harga eigen karakteristik adalah sebagai berikut :
J X X
(19b)
J merupakan matrik persegi N N dan merupakan harga eigen
dimana
karakteristik yang dapat ditentukan dengan menyelesaikan persamaan karakteristik berikut (lihat persamaan 10 hal. 8) 11,12
det J I 0 ,
dimana
I adalah
matriks identitas
3 x 2 2ax 1 1 0 0 det 2dx 1 0 0 0 0 0 0 0 b det
3 x 2 2ax
1
2dx b
1 0
1 0 0 (20)
sesuai dengan aturan Cramer maka didapatkan persamaan polinomial untuk sebagai berikut :
20
(3 x 2 2ax ) ( 1 ) ( ) 0 (1)(1 )( b) ( )(1)(2dx)
(21)
untuk I 1.78 , a 3 , b 4, d 5 ,
0.001 , dan xe
1 2
5 1 ,
(b)
Maka didapatkan akar-akar dari harga eigen karakteristik sebagai berikut:
Gambar 9. Analisa bifurkasi pada persamaan H-R
1 -12.7193 , 2 0.0166 + 0.0018i , 3 0.0166 - 0.0018i.
Pada Gambar 9a merupakan pola potensial aku si yang bentuk dari pola bifurkasi pada trayektori Gambar 9b dimana pada Gambar 9 pola trayektori membentuk suatu kecenderungan bifurkasi yang terjadi dari sistem dan membentuk bifurkasi Hopf dan bifurkasi Sadle-Node. Bifurkasi Hopf pada dasarnya melibatkan trayektori yang bersifat periodik dimana terjadi perubahan jenis titik kritis dari titik fokus dengan atraktor positif menjadi atraktor negatif disertai dengan kemunculan limit cycle. Limit cycle merupakan sebuah trayektori berbentuk lingkaran yang bersifat periodik yang muncul akibat perubahan kestabilan fokus. Sedangkan bifurkasi Sadle-Node yang mungkin terjadi akibat munculnya dua atau lebih dari titik kritis. Dalam kasus ini jumlah Limit cycle akan menentukan pola spiking yang terjadi didalam sistem persamaan H-R.
Harga-harga eigen karakteristik diatas menentukan jenis-jenis bifurkasi pada dinamika sistem penjalaran impuls pada sel saraf. Dengan menggunakan Tabel 2 (hal. 8) harga eigen yang dimiliki sistem saat I 1.78 adalah terdiri dari bilangan rill dan kompleks sehingga titik fokus terbentuk pada ruang fasa. Pada perubahan input berupa I akan merubah pola trayektori yang terjadi pada sistem pola tersebut dijelaskan pada Gambar 9a dan Gambar 9b.
(a)
Tabel 4 menunjukkan beberapa titik yang mengidentifikasi bifurkasi yang terjadi pada dinamika sistem penjalaran impuls di dalam sel saraf berdasarkan model H-R. Titik-titik tersebut dihasilkan dari perubahan input berupa I arus searah yang diubah untuk menemukan beberapa titik kritis pada sistem dan menemukan karakteristik dari trayektori sistem. Pola bifurkasi berperan penting dalam menemukan sejumlah fenomena dalam penjalaran impuls saraf atau pembentukan pola potensial aksi selama menjalar sepanjang sel saraf.
21
Tabel 4. Bifurkasi yang terjadi pada dinamika sistem dengan variasi I input
I
x0 y0 z 0
1 2 3
Jenis bifurkasi
0.81
1.45 9.57 0.66)
16.03 0.02 0.01
SadleNode
1.78
1.21 6.35 1.62
12.72 0.02 0.002i 0.02 0.002i
Hopf
3.32
0.69 1.42 3.68
6.82 0.19 0.002
SadleNode
y
(a)
Pada titik tertentu dalam trayektori membentuk suatu pola bifurkasi seperti tertera pada Gambar 10a dan Gambar 10c adalah bifurkasi Sadle-Node terbentuk dengan melihat harga eigen karakteristik bernilai rill dan lebih dari satu. Sementara itu bifurkasi Hopf Gambat 10b terbentuk dari adanya nilai eigen karakteristik yang bernilai imajiner. Letak titik kesetimbangan pada sistem selalu berubah dengan perubahan nilai input I berupa arus searah yang akan menjadikan sistem membentuk suatu pola potensial aksi yang menjalar sepanjang sel saraf. (b)
22
Limit cycle
(c) Gambar 10. Analisis bifurkasi dengan variasi I input. (a) Bifurkasi SadelNode I 0.81 . (b) Bifurkasi Hopf I 1.78 . (c) Bifurkasi Sadel-Node
I 5.32
Analisa dinamika sistem pada penjalaranan impuls di dalam sel saraf akan melihat lebih luas dari yang sebenarnya terjadi pada persamaan. Dengan adanya titik-titik kesimbangan dan jenis bifurkasi yang menyertainya sistem tersebut berubah terhadap variasi kemudian perubahan I inputan, trayektori yang menunjukkan bagaimana persamaan membangun suatu dinamika sistem pada sel saraf.
4.3 GUI untuk Dinamika Sistem H-R dengan MATLAB Graphical User Interface atau yang disingkat GUI merupakan suatu sarana singkat dan populer untuk menampilkan aplikasi perangkat lunak sehingga lebih mudah digunakan. Didalam pembuatan GUI para programmer (pembuat perangkat lunak) haruslah
memperhitungkan segi kenyamanan pengguna untuk mengakses perangkat lunak tersebut. Sehingga apa yang dimaksud dari aplikasi yang dibuat akan mudah dimengerti oleh pengguna. Dinamika sistem penjalaran impuls pada sel saraf dengan menggunakan model HR memiliki banyak aspek dan variabel yang ditinjau. Pada sub-bab sebelumnya telah dibahas mengenai solusi numerik yang menampilkan bentuk fisis dari persamaan H-R dan kemudian telah dibahas mengenai dinamika sistem serta bifurkasinya. Namun kesemua hal tersebut hanya dapat diakses dan dipahami oleh pihak-pihak tertentu. GUI kali ini akan menampilkan suatu aplikasi tentang simulasi penjalaran impuls pada sel saraf dengan menggunakan model HR. Pada Gambar 11 diperlihatkan fisik dari aplikasi yang dimaksud. Pada aplikasi tersebut terdapat variabel, konstanta, dan inputan yang menyusun persamaan H-R, serta terdapat pula keluaran apabila program dijalankan dengan menekan tombol ‘PROSES’ dan hasil akan ditampilkan pada axis pertama (lihat Gambar 11). Tombol tersbut merupakan pemicu program yang akan menjalankan serangkaian proses simulasi dengan metode RKF45 yang telah disebutkan sebelumnya dari model H-R. Bila telah menampakkan hasil pada grafik keluaran yang pertama maka dapat menganalisa tindajaun berikutnya dengan menekan tombol ‘PLOT’. Penggunaan tombol ‘PLOT’ adalah untuk mengakses kotak edit dibawah label ‘ZOOM’ sehingga semua axis akan terisi dengan keluaran yang diinginkan serta dapat dianalisis dengan melihat label dari setiap axis (lihat Gambar 12). Sedangkan tomnol ‘RESET’ adalah tombol untuk mengembalikan semua kondisi aplikasi setelah melakukan perubahan terhadapat variabel, konstanta, dan inputan dari persamaan H-R.
23
Gambar 11. GUI simulasi H-R model. Sejumlah kontanta dan input dapat diubah untuk melakukan analisis dari penjalaran impuls pada sel saraf. Keluaran pertama setelah menekan tombol ‘PROSES’.
Gambar 12. ‘Penggunaan tombol ‘PLOT’. Untuk mendapatkan grafik yang lebih dekat sehingga memudahkan untuk dianalisis.
24
Pada sub bab sebelumnya dibahas mengenai konstanta xe yang merupakan titik keseimbangan saat sistem jika tidak diberi arus atau impuls. Gambar 8 pada GUI disajikan saat I 0 maka akan terbentuk titik keseimbangan pada x0 xe hal tersebut berarti pada
kondisi tidak diberi arus sistem akan berada dalam keseimbangan titik-titik tersebut dan inilah merupakan fisis dari sifat sel saraf yang selalu menuju pada keadaan seimbang sesaat setelah diberi rangsangan berupa arus atau impuls.
Gambar 13. ‘Sifat dari konstanta xe sebagai titik kesembangan bagi x0 pada saat I 0 pada sistem’. Nilai tersebut menjelaskan sifat sel saraf yang selalu menuju ke keadaan semula sesaat setelah diberi rangsangan.
BAB V KESIMPULAN DAN SARAN 5. 1 Kesimpulan Simulasi dinamika sel saraf menggunakan model Hindmarsh-Rose merupakan suatu simulasi yang dapat dilakukan untuk menjelaskan tentang dinamika penjalaran yang terjadi didalam sel saraf. Model ini dapat menunjukkan fenomena-fenomena seperti spiking-bursting dan chaotic bursting dengan mengubah input berupa arus searah yang merupakan fenomena pengembangan dari model penjalaran impuls didalam sel saraf. Perubahan
input berupa arus searah akan mengakibatkan transisi fenomena dari spike menjadi bursting dan memunculkan fenomena chotic bursting sepanjang penjalaran impuls didalam sel saraf. Fenomena chotic bursting memungkinkan terbentuknya fase peningkatan jumlah spike untuk kenaikan 0.05A dari input kedalam sistem. Peningkatan tersebut terjadi secara linier pada range input I 3.27A dan menampakkan fenomena menyimpang pada range diatas input tersebut oleh keterbatasan dari persamaan dan perilaku dari sel saraf yang sebenarnya. Metode RKF 45 merupakan metode numerik untuk memecahkan persamaan
25
differensial berorde yang terdapat dalam sistem persamaan Hindmarsh-Rose. Analisa dinamika sistem melibatkan titik-titik keseimbangan dan bifurkasi lokal pada persamaan non-linier yang terdapat didalam sistem tersebut. Dengan menemukan titik-titik tersebut akan lebih mempertajam penjelasan mengenai dinamika sistem yang dibentuk dari persamaan tersebut. Simulasi berbasis GUI yang diimplementasikan pada perangkat lunak MATLAB merupakan suatu sajian yang sederhana dan user friendly untuk melakukan analisa terhadap dinamika penjalaran impuls pada sel saraf. Tampilan berbasis GUI akan mengharmoniskan hubunngan inputoutput dari suatu masalah sehingga lebih terintegrasi dan layak untuk di publikasikan. Dari keunggulan tersebut hasil simulasi ini dapat digunakan oleh lebih banyak pengguna tanpa harus memahami metode numerik lebih yang lebih sulit.
5.2 Saran Penelitian ini merupakan penelitian dasar dengan hanya melakukan tinjuan analitik dan membuat simulasi sederhana, sehingga untuk mengembangkan kearah sistem dinamika yang lebih kompleks seperti membuat pengkopelan terhadap beberapa sel saraf. Spesifikasi dari perangkat keras dan faktor eksternal dari simulasi ini bisa saja terjadi diluar kendali dari peneliti. Untuk simulasi berbasis GUI bisa dibuat kembali dengan lebih mempertimbangkan taraf kemudahan bagi pengguna seperti tampilan output berupa indeks dan legend dari pada axis hingga peletakan tombol, slide bar, dan edit yang lebih tepat. Bagian analisa pada persamaan ini perlu ditinjau lagi dengan mengganti beberapa konstanta sehingga dapat terlihat bagaimana perilaku persamaan pada kondisi berbeda. Analisa dinamika sistem yang terdapat pada persamaan ini bisa ditelaah lebih dalam sebagai acuan untuk mememukan ataupun menyempurnakan model-model dinamika yang terdapat pada bidang biofisika.
26
DAFTAR PUSTAKA 1
2
3
4
5
6
Tipler, P. A., 2003, FISIKA Untuk Sains dan Teknik Jilid 2, Penerbit Erlangga, Jakarta, , hal 163-167. Neefs, P. J, 2008, On Time-Delayed Couple Hindmarsh-Rose Neurons: Stability of Equilibria and Eindhoven Synchronosation, University of Technology, DCT, hal 1-12. Hindmarsh J. L, 1982, Rose R. M, A Model of The Nerve Impulse Using Two First-Order Differential Equations, Cardiff CF1 1XL, UK, 1982 Nature Vol.296 hal 162-164. Fitzhugh, R., 1961, Impulses and Physiological States in Theoretical Models of Nerve Membrane, National Insitutes of Health, Bathesda, Biophysical Journal Vol. 1. Corson N. dan Aziz-Alaoui M. A., 2006, Dynamics and complexity of Hinmarsh-Rose neural systems, Jurnal, Le Harve. Adi. Neuron. Januari 2008. Web. 16 November 2010.
.
7
Izhikevich E. M, 2007, Dynamical Systems in Neuroscience The Geometry of Excitability and Bursting, The MIT Press, London, hal 25-27.
8
Yeargers, K. Edward,dkk., 1996, An Introduction to the mathematics of Biology, Birkhauser, Boston, hal 15
9
Dinamika Yuarnita, D., 2009, Iimpuls Pada Satu Sel Saraf dengan Akupuntur Sebagai Stimulus, Institut Teknologi Bandung, Skripsi, hal 14.
10 (Anonim). Dinamika Sistem. Web. 30 November 2010. <www.wikipedia.org/wiki/Dinamik a_Sistem.html> 11 Saputra, A., Solusi Periodik Eliptik Jacobi Persamaan Modus Tergandeng Sistem Kisi Bragg,Departemen Fisika FMIPAIPB,Skripsi,2010, hal 3-6. 12 Alatas, H., 2010, Modul Fisika Non-liniear , Departemen Fisika FMPA-IPB, Bogor, Bab III Dinamika Sistem. 13 Hillborn, R. C., Chaos and Oxford Nonlinear Dynamics, University Press, USA, 14 Steven H. Strogatz, 1994, Nonlinear Dynamics and Chaos, Perseus Books Publishing, L.L.C.,New York,. 15 Mathews , J. H., Fink, K. K., 2003, Numerical Methods Using Matlab, 3th Edition, Prentice-Hall Inc. Upper Saddle River, New Jersey, USA,hal 497-499. 16 Supriyanto, 2007, Komputasi untuk Sains dan Teknik, Universitas Indonesia, Jakarta , Bab 10 Persamaan Differensial. 17 Bismo, S, 2009, Seri Mata Kuliah : Permodelan dan Matematika Terapan, Universitas Indonesia, Jakarta, Modul 5. 18 Away, Gunaidi Abdia, 2010. The Shortcut of MATLAB Bandung: Programming. Informatika.
27
LAMPIRAN
28
Lampiran 1. Diagram Alir Penelitian
Mulai
Penelusuran literatur
Pembuatan program dengan metode numerik RKF45 dari persamaan Hindmarsh-Rose (perangkat lunak MATLAB)
Menemukan solusi numerik persamaan Hindmarsh-Rose
Pembuatan Graphical User Interface (GUI)
Analisa simulasi dinamika sistem Variasi sejumlah inputan untuk mendapatkan fenomena struktur dari dinamika sistem sel saraf Menemukan titik kritis dan jenis bifurkasi lokal Menganalisa tingkah laku dinamika sistem pada sel saraf
Membandingkan dengan literatur
Pembahasan dan penyusunan hasil peneliatian
Selesai
29
Lampiran 2. Persamaan Hodgkin-Huxley untuk Potensial Aksi pada Satu Sel Saraf
3 4 dV C Na m h( E Na V ) C K n ( E K V ) dt Cleak (Vrest V ) I inj (t )
Potensial membran
:
Sodium I Na Fast
:
dm m (V )(1 m) m (V )m dt
Sodium I Na Slow
:
dh h (V )(1 h) h (V )h dt
Potasium I Na Slow
:
dn n (V )(1 n) n (V )n dt
n (V )
0.0110 V ; n 0.125 exp V / 80 10 V exp 1 10
h (V ) 0.07 exp V / 20 ; h (V )
m (V )
1 30 V exp 1 10
0.125 V ; m (V ) 4 exp V / 18 25 V exp 1 10
30
Lampiran 3. Program Matlab untuk Metode RKF 45 dari persamaan Hindmarsh-Rose %================================================================= %SIMULASI DINAMIKA SISTEM SEL SARAF MENGGUNAKAN MODEL HINDMARSHROSE %================================================================= %Persamaan H-R % x(.) = y+a*x^2-x^3-z+I % y(.) = 1-d*x^2-y % z(.) = m*(b*(x-xe)-z) %----------------------------------------------------------------clear all clc % Penyediaan matriks x=[]; y=[]; z=[]; t=[]; % Konstanta model Hindmarsh-Rose a=3; b=4; d=5; xe=-1/2*(1+sqrt(5)); m=0.001; % Masukan arus (rangsangan berupa arus listrik) I=3.6; % Inisial kondisi xi=0; yi=0; zi=0; % Waktu dan Jumlah iterasi (Pencacahan) ti=0; tf=5000; n=50000; h=(tf-ti)/n; % Pengisian inisial matriks x(1)=xi; y(1)=yi; z(1)=zi; t(1)=ti; % Iterasi Model H-R dengan metode RKF45 for k=1:n kx1=h*(y(k)+a*x(k)^2-x(k)^3-z(k)+I); ky1=h*(1-d*(x(k))^2-(y(k))); kz1=h*(m*(b*(x(k)-xe)-z(k))); kx2=h*((y(k)+1/4*ky1)+a*(x(k)+1/4*kx1)^2-(x(k)+1/4*kx1)^3(z(k)+1/4*kz1)+I); ky2=h*(1-d*(x(k)+1/4*kx1)^2-(y(k)+1/4*ky1)); kz2=h*(m*(b*((x(k)+1/4*kx1)-xe)-(z(k)+1/4*kz1))); kx3=h*((y(k)+3/32*ky1+9/32*ky2)+a*(x(k)+3/32*kx1+9/32*kx2)^2(x(k)+3/32*kx1+9/32*kx2)^3-(z(k)+3/32*kz1+9/32*kz2)+I); ky3=h*(1-d*(x(k)+3/32*kx1+9/32*kx2)^2-(y(k)+3/32*ky1+9/32*ky2)); kz3=h*(m*(b*((x(k)+3/32*kx1+9/32*kx2)-xe)(z(k)+3/32*kz1+9/32*kz2))); kx4=h*((y(k)+1932/2197*ky17200/2197*ky2+7296/2197*ky3)+a*(x(k)+1932/2197*kx17200/2197*kx2+7296/2197*kx3)^2-(x(k)+1932/2197*kx17200/2197*kx2+7296/2197*kx3)^3-(z(k)+1932/2197*kz17200/2197*kz2+7296/2197*kz3)+I); ky4=h*(1-d*(x(k)+1932/2197*kx1-7200/2197*kx2+7296/2197*kx3)^2(y(k)+1932/2197*ky1-7200/2197*ky2+7296/2197*ky3)); kz4=h*(m*(b*((x(k)+1932/2197*kx1-7200/2197*kx2+7296/2197*kx3)-xe)(z(k)+1932/2197*kz1-7200/2197*kz2+7296/2197*kz3))); kx5=h*((y(k)+439/216*ky1-8*ky2+3680/513*ky3845/4104*ky4)+a*(x(k)+439/216*kx1-8*kx2+3680/513*kx3845/4104*kx4)^2-(x(k)+439/216*kx1-8*kx2+3680/513*kx3845/4104*kx4)^3-(z(k)+439/216*kz1-8*kz2+3680/513*kz3845/4104*kz4)+I); ky5=h*(1-d*(x(k)+439/216*kx1-8*kx2+3680/513*kx3-845/4104*kx4)^2(y(k)+439/216*ky1-8*ky2+3680/513*ky3-845/4104*ky4)); kz5=h*(m*(b*((x(k)+439/216*kx1-8*kx2+3680/513*kx3-845/4104*kx4)xe)-(z(k)+439/216*kz1-8*kz2+3680/513*kz3-845/4104*kz4)));
31
kx6=h*((y(k)-8/27*ky1+2*ky2-3544/2565*ky3+1859/4104*ky411/40*ky5)+a*(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^2-(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^3-(z(k)-8/27*kz1+2*kz2-3544/2565*kz3+1859/4104*kz411/40*kz5)+I); ky6=h*(1-d*(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^2-(y(k)-8/27*ky1+2*ky2-3544/2565*ky3+1859/4104*ky411/40*ky5)); kz6=h*(m*(b*((x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)-xe)-(z(k)-8/27*kz1+2*kz2-3544/2565*kz3+1859/4104*kz411/40*kz5))); t(k+1)=t(k)+h; x(k+1)=x(k)+16/135*kx1+6656/12825*kx3+28561/56430*kx49/50*kx5+2/55*kx6; y(k+1)=y(k)+16/135*ky1+6656/12825*ky3+28561/56430*ky49/50*ky5+2/55*ky6; z(k+1)=z(k)+16/135*kz1+6656/12825*kz3+28561/56430*kz49/50*kz5+2/55*kz6; end % Plot data plot(t,x) xlabel('t (detik)') ylabel('x (satuan potensial)') figure plot3(x,y,z) xlabel('x') ylabel('y') zlabel('z') % Pencarian titik keseimbangan sistem x1=[1 (d-a) b (-b*xe-I-1)]; x0=roots(x1) for i=1:3 y0(i)=1-d*x0(i)^2; z0(i)=(b*(x0(i)-xe)); end y0=y0' z0=z0' % Pencarian karaktekterisasi titik keseimbangan sistem jacob1=[-3*x0(1)^2+2*a*x0(1) 1 -1;-2*d*x0(1) -1 0;m*b 0 -m] jacob2=[-3*x0(2)^2+2*a*x0(2) 1 -1;-2*d*x0(2) -1 0;m*b 0 -m] jacob3=[-3*x0(3)^2+2*a*x0(3) 1 -1;-2*d*x0(3) -1 0;m*b 0 -m] eigent1=eig(jacob1) eigent2=eig(jacob2) eigent3=eig(jacob3)
32
Lampiran 4. Program Matlab untuk GUI Fisik HR_Simulator.m clear all; clc; ruang=figure('units','pixel','position',[15 25 1250 700],'color',[.8 .8 .8],'menubar','figure','resize','on','name','H-R Simulator','NumberTitle','off'); judul=uicontrol('parent',ruang,'units','pixel','position',[50 650 900 35],'style','text','string','HINDMARSH-ROSE NEURON MODEL SIMULATOR','backgroundcolor',[0.8 0.8 0.8],'fontname','AddElectricCity','fontsize',24,'fontweight','bold ','foregroundcolor',[0 0 0],'horizontalalignment','center'); grafik1=axes('parent',ruang,'units','pixel','position',[220 400 350 165],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); grafik2=axes('parent',ruang,'units','pixel','position',[220 175 350 165],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); grafik3=axes('parent',ruang,'units','pixel','position',[630 470 165 155],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); grafik4=axes('parent',ruang,'units','pixel','position',[630 280 165 155],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); grafik5=axes('parent',ruang,'units','pixel','position',[630 90 165 155],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); grafik6=axes('parent',ruang,'units','pixel','position',[865 180 350 350],'xgrid','on','ygrid','on','xcolor',[.4 0 .15],'ycolor',[.4 0 .15],'fontsize',8,'color',[1 1 1]); tomproses=uicontrol('parent',ruang,'units','points','position',[23 0 30 60 25],'style','pushbutton','callback','HR_Rk45','string','PROSES','f ontname','Tahoma','fontweight','bold','fontsize',10); tomplot=uicontrol('parent',ruang,'units','points','position',[300 30 60 25],'style','pushbutton','callback','HR_Plot','string','PLOT','fon tname','Tahoma','fontweight','bold','fontsize',10); tomreset=uicontrol('parent',ruang,'units','points','position',[370 30 60 25],'style','pushbutton','callback','HR_Reset','string','RESET','f ontname','Tahoma','fontweight','bold','fontsize',10); label1=uicontrol('parent',ruang,'units','pixel','position',[20 600 30 20],'style','text','string','a=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label2=uicontrol('parent',ruang,'units','pixel','position',[20 570 30 20],'style','text','string','b=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label3=uicontrol('parent',ruang,'units','pixel','position',[95 600 30 20],'style','text','string','d=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label4=uicontrol('parent',ruang,'units','pixel','position',[95 570 30 20],'style','text','string','I=','backgroundcolor',[.8 .8
33
.8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label5=uicontrol('parent',ruang,'units','pixel','position',[20 530 30 20],'style','text','string','m=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label6=uicontrol('parent',ruang,'units','pixel','position',[15 500 35 20],'style','text','string','xe=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label7=uicontrol('parent',ruang,'units','pixel','position',[15 460 35 20],'style','text','string','ti=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label8=uicontrol('parent',ruang,'units','pixel','position',[15 430 35 20],'style','text','string','tf=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label9=uicontrol('parent',ruang,'units','pixel','position',[15 50 35 20],'style','text','string','xi=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label10=uicontrol('parent',ruang,'units','pixel','position',[100 50 35 20],'style','text','string','yi=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label11=uicontrol('parent',ruang,'units','pixel','position',[180 50 35 20],'style','text','string','zi=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label12=uicontrol('parent',ruang,'units','pixel','position',[15 390 35 20],'style','text','string','xo=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label13=uicontrol('parent',ruang,'units','pixel','position',[15 360 35 20],'style','text','string','yo=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label14=uicontrol('parent',ruang,'units','pixel','position',[15 330 35 20],'style','text','string','zo=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label15=uicontrol('parent',ruang,'units','pixel','position',[15 290 35 20],'style','text','string','e1=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label16=uicontrol('parent',ruang,'units','pixel','position',[15 265 35 20],'style','text','string','e2=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label17=uicontrol('parent',ruang,'units','pixel','position',[15 235 35 20],'style','text','string','e3=','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right'); label18=uicontrol('parent',ruang,'units','pixel','position',[80 165 45 40],'style','text','string','ZOOM','backgroundcolor',[.8 .8 .8],'fontname','AddElectricCity','fontsize',12,'fontweight','Norma l','foregroundcolor',[0 0 0],'horizontalalignment','right');
34
isi1=uicontrol('parent',ruang,'units','pixel','position',[55 600 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','3','fontname','AddElectricCity','fonts ize',8); isi2=uicontrol('parent',ruang,'units','pixel','position',[55 570 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','4','fontname','AddElectricCity','fonts ize',8); isi3=uicontrol('parent',ruang,'units','pixel','position',[130 600 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','5','fontname','AddElectricCity','fonts ize',8); isi4=uicontrol('parent',ruang,'units','pixel','position',[130 570 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','3.25','fontname','AddElectricCity','fo ntsize',8); isi5=uicontrol('parent',ruang,'units','pixel','position',[55 530 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0.001','fontname','AddElectricCity','f ontsize',8); isi6=uicontrol('parent',ruang,'units','pixel','position',[55 500 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','1.6180','fontname','AddElectricCity','fontsize',8); isi7=uicontrol('parent',ruang,'units','pixel','position',[55 460 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi8=uicontrol('parent',ruang,'units','pixel','position',[55 430 40 20],'backgroundcolor',[1 1 1],'style','Edit','string','3000','fontname','AddElectricCity','fo ntsize',8); isi9=uicontrol('parent',ruang,'units','pixel','position',[100 445 50 20],'backgroundcolor',[1 1 1],'style','Edit','string','50000','fontname','AddElectricCity','f ontsize',8); isi10=uicontrol('parent',ruang,'units','pixel','position',[55 50 40 20],'backgroundcolor',[.8 .8 .8],'style','Edit','string','0','fontname','AddElectricCity','font size',8); isi11=uicontrol('parent',ruang,'units','pixel','position',[140 50 40 20],'backgroundcolor',[.8 .8 .8],'style','Edit','string','0','fontname','AddElectricCity','font size',8); isi12=uicontrol('parent',ruang,'units','pixel','position',[225 50 40 20],'backgroundcolor',[.8 .8 .8],'style','Edit','string','0','fontname','AddElectricCity','font size',8); isi13=uicontrol('parent',ruang,'units','pixel','position',[55 390 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi14=uicontrol('parent',ruang,'units','pixel','position',[55 360 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi15=uicontrol('parent',ruang,'units','pixel','position',[55 330 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8);
35
isi16=uicontrol('parent',ruang,'units','pixel','position',[55 295 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi17=uicontrol('parent',ruang,'units','pixel','position',[55 265 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi18=uicontrol('parent',ruang,'units','pixel','position',[55 235 80 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi19=uicontrol('parent',ruang,'units','pixel','position',[55 165 50 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi20=uicontrol('parent',ruang,'units','pixel','position',[55 135 50 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi21=uicontrol('parent',ruang,'units','pixel','position',[110 165 50 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8); isi22=uicontrol('parent',ruang,'units','pixel','position',[110 135 50 20],'backgroundcolor',[1 1 1],'style','Edit','string','0','fontname','AddElectricCity','fonts ize',8);
36
Lampiran 5. Program Matlab untuk GUI Isi Program HR_Rk45.m %-------------------------------------------------------------%dx/dt=y+a*x^2-x^3-z+I %dy/dt=1-d*x^2-y %dz/dt=m*(b*(x-xe)-z) %-------------------------------------------------------------a=str2num(get(isi1,'string')); b=str2num(get(isi2,'string')); d=str2num(get(isi3,'string')); I=str2num(get(isi4,'string')); m=str2num(get(isi5,'string')); xe=str2num(get(isi6,'string')); ti=str2num(get(isi7,'string')); tf=str2num(get(isi8,'string')); n=str2num(get(isi9,'string')); xi=str2num(get(isi10,'string')); yi=str2num(get(isi11,'string')); zi=str2num(get(isi12,'string')); x=[]; y=[]; z=[]; t=[]; h=(tf-ti)/n; x(1)=xi; y(1)=yi; z(1)=zi; t(1)=ti; for k=1:n kx1=h*(y(k)+a*x(k)^2-x(k)^3-z(k)+I); ky1=h*(1-d*(x(k))^2-(y(k))); kz1=h*(m*(b*(x(k)-xe)-z(k))); kx2=h*((y(k)+1/4*ky1)+a*(x(k)+1/4*kx1)^2-(x(k)+1/4*kx1)^3(z(k)+1/4*kz1)+I); ky2=h*(1-d*(x(k)+1/4*kx1)^2-(y(k)+1/4*ky1)); kz2=h*(m*(b*((x(k)+1/4*kx1)-xe)-(z(k)+1/4*kz1))); kx3=h*((y(k)+3/32*ky1+9/32*ky2)+a*(x(k)+3/32*kx1+9/32*kx2)^2(x(k)+3/32*kx1+9/32*kx2)^3-(z(k)+3/32*kz1+9/32*kz2)+I); ky3=h*(1-d*(x(k)+3/32*kx1+9/32*kx2)^2-(y(k)+3/32*ky1+9/32*ky2)); kz3=h*(m*(b*((x(k)+3/32*kx1+9/32*kx2)-xe)(z(k)+3/32*kz1+9/32*kz2))); kx4=h*((y(k)+1932/2197*ky17200/2197*ky2+7296/2197*ky3)+a*(x(k)+1932/2197*kx17200/2197*kx2+7296/2197*kx3)^2-(x(k)+1932/2197*kx17200/2197*kx2+7296/2197*kx3)^3-(z(k)+1932/2197*kz17200/2197*kz2+7296/2197*kz3)+I); ky4=h*(1-d*(x(k)+1932/2197*kx1-7200/2197*kx2+7296/2197*kx3)^2(y(k)+1932/2197*ky1-7200/2197*ky2+7296/2197*ky3)); kz4=h*(m*(b*((x(k)+1932/2197*kx1-7200/2197*kx2+7296/2197*kx3)-xe)(z(k)+1932/2197*kz1-7200/2197*kz2+7296/2197*kz3))); kx5=h*((y(k)+439/216*ky1-8*ky2+3680/513*ky3845/4104*ky4)+a*(x(k)+439/216*kx1-8*kx2+3680/513*kx3845/4104*kx4)^2-(x(k)+439/216*kx1-8*kx2+3680/513*kx3-
37
845/4104*kx4)^3-(z(k)+439/216*kz1-8*kz2+3680/513*kz3845/4104*kz4)+I); ky5=h*(1-d*(x(k)+439/216*kx1-8*kx2+3680/513*kx3-845/4104*kx4)^2(y(k)+439/216*ky1-8*ky2+3680/513*ky3-845/4104*ky4)); kz5=h*(m*(b*((x(k)+439/216*kx1-8*kx2+3680/513*kx3-845/4104*kx4)xe)-(z(k)+439/216*kz1-8*kz2+3680/513*kz3-845/4104*kz4))); kx6=h*((y(k)-8/27*ky1+2*ky2-3544/2565*ky3+1859/4104*ky411/40*ky5)+a*(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^2-(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^3-(z(k)-8/27*kz1+2*kz2-3544/2565*kz3+1859/4104*kz411/40*kz5)+I); ky6=h*(1-d*(x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)^2-(y(k)-8/27*ky1+2*ky2-3544/2565*ky3+1859/4104*ky411/40*ky5)); kz6=h*(m*(b*((x(k)-8/27*kx1+2*kx2-3544/2565*kx3+1859/4104*kx411/40*kx5)-xe)-(z(k)-8/27*kz1+2*kz2-3544/2565*kz3+1859/4104*kz411/40*kz5))); t(k+1)=t(k)+h; x(k+1)=x(k)+16/135*kx1+6656/12825*kx3+28561/56430*kx49/50*kx5+2/55*kx6; y(k+1)=y(k)+16/135*ky1+6656/12825*ky3+28561/56430*ky49/50*ky5+2/55*ky6; z(k+1)=z(k)+16/135*kz1+6656/12825*kz3+28561/56430*kz49/50*kz5+2/55*kz6; end set(ruang,'currentaxes',grafik1); plot(t,x) xlabel('t') ylabel('x') set(isi19,'string',num2str(ti+1)); set(isi21,'string',num2str(n)); set(isi20,'string',num2str(ti+1)); set(isi22,'string',num2str(n));
x1=[1 (d-a) b (-b*xe-I-1)]; x0=roots(x1); set(isi13,'string',num2str(x0(3))); y0=1-d*x0(3)^2; z0=(b*(x0(3)-xe)); set(isi14,'string',num2str(y0)); set(isi15,'string',num2str(z0)); jacob=[-3*x0(3)^2+2*a*x0(3) 1 -1;-2*d*x0(3) -1 0;m*b 0 -m]; eigent=eig(jacob); set(isi16,'string',num2str(eigent(1))); set(isi17,'string',num2str(eigent(2))); set(isi18,'string',num2str(eigent(3)));
38
HR_Plot.m ai=str2num(get(isi19,'string')); af=str2num(get(isi21,'string')); bi=str2num(get(isi20,'string')); bf=str2num(get(isi22,'string'));
set(ruang,'currentaxes',grafik1); plot(t(1,(ai:af)),x(1,(ai:af))); xlabel('t'); ylabel('x'); set(ruang,'currentaxes',grafik2); plot(t(1,(bi:bf)),x(1,(bi:bf))); xlabel('t'); ylabel('x'); set(ruang,'currentaxes',grafik3); plot(x(1,(bi:bf)),y(1,(bi:bf))); xlabel('x'); ylabel('y'); set(ruang,'currentaxes',grafik4); plot(x(1,(bi:bf)),z(1,(bi:bf))); xlabel('x'); ylabel('z'); set(ruang,'currentaxes',grafik5); plot(y(1,(bi:bf)),z(1,(bi:bf))); xlabel('y'); ylabel('z'); set(ruang,'currentaxes',grafik6); plot3(x(1,(bi:bf)),y(1,(bi:bf)),z(1,(bi:bf))); xlabel('x'); ylabel('y'); zlabel('z');
HR_Reset.m set(isi1,'string',num2str(3)); set(isi2,'string',num2str(4)); set(isi3,'string',num2str(5)); set(isi4,'string',num2str(3.25)); set(isi5,'string',num2str(0.001)); set(isi6,'string',num2str(-1.6180)); set(isi7,'string',num2str(0)); set(isi8,'string',num2str(3000)); set(isi9,'string',num2str(50000)); set(isi10,'string',num2str(0)); set(isi11,'string',num2str(0)); set(isi12,'string',num2str(0)); set(isi13,'string',num2str(0)); set(isi14,'string',num2str(0)); set(isi15,'string',num2str(0)); set(isi16,'string',num2str(0)); set(isi17,'string',num2str(0)); set(isi18,'string',num2str(0)); set(isi19,'string',num2str(0)); set(isi20,'string',num2str(0)); set(isi21,'string',num2str(0)); set(isi22,'string',num2str(0));
