Mètodes de control discrets per alternants cardíacs

Títol: Mètodes de control discrets per alternants cardíacs Autora: Nora Wieczorek i Masdeu Director: Blas Echebarria Departament: Departament de Física Convocatòria: 2017-2018 Grau en Matemàtiques Universitat Polit` ecnica de Catalunya Facultat de Matem` atiques i Estad´ ıstica Treball de Final de Grau M`etodes de control discrets per alternans card´ıacs Nora Wieczorek i Masdeu Tutor del treball Blas Echebarria 4 de setembre de 2018 Agra¨ıments Vull agrair al Blas la seva dedicaci´o de temps i el suport que m’ha donat durant aquest treball. 3 Abstract Els alternans s´on variacions peri`odiques en el potencial d’acci´o de les c`el·lules card´ıaques. Aquests poden provocar el pas del ritme normal del cor a la taquic`ardia e incl`os a la fibril·laci´o, amb la p`erdua de la capacitat de bombeig del cor, que sovint resulta en mort card´ıaca sobtada. En aquest treball, tractem els alternans tant a nivell unicel·lular com des d’un teixit unidimensional d’una certa longitud, el qual descriu el comportament que segueix un seguit de c`el·lules tenint en compte els seus enlla¸cos. En primer lloc, estudiem el perqu`e de la aparici´o d’alternans i plantejem dos m`etodes de control discrets en els que la variable ser`a el per´ıode de batec del cor. Hem vist que la correcci´o ´es possible per tot per´ıode si els m`etodes segueixen una s`erie de condicions. En el teixit card´ıac, apliquem les correccions del cas unicel·lular i observem l’efecte que tenen els m`etodes al extendre’ls a teixit. En aquest, hem vist que la variable longitud tamb´e repercuteix en l’efic`acia dels controls que funcionen de manera similar: deixen de controlar un cop arribada una longitud m`axima de teixit. Paraules Clau: Biologia matem`atica, Biof´ısica, Din`amica card´ıaca, Alternans card´ıacs, M`etodes de control discrets, Estabilitat de sistemes discrets, M`etodes num`erics. 5 ´ Index 1 Introducci´o 10 1.1 Fisiologia a nivell cel·lular............................... 11 1.2 Propagaci´o de l’ona el`ectrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 Alternans, el fenomen que volem corregir . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Objectiusdeltreball .................................. 15 2 Model matem`atic 16 2.1 Model d’una sola c`el·lula (o zero dimensional) . . . . . . . . . . . . . . . . . . . . 16 2.1.1 Sistemaexcitable................................ 21 2.1.2 Resoluci´o num`erica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 Model del teixit card´ıac (unidimensional) . . . . . . . . . . . . . . . . . . . . . . 23 2.2.1 Resoluci´onum`erica............................... 24 3 La corba de restituci´o i la seva estabilitat 27 3.1 Conceptes te`orics de sistemes din`amics . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Plantejament del sistema: la corba de restituci´o . . . . . . . . . . . . . . . . . . . 28 3.2.1 C`alcul num`eric de punts de la corba de restituci´o . . . . . . . . . . . . . 29 3.2.2 C`alcul anal´ıtic de la corba de restituci´o . . . . . . . . . . . . . . . . . . . 30 3.3 Estabilitatdelsistema................................. 32 3.4 Resultatsnum`erics................................... 34 3.4.1 Resultats num`erics pel model zero dimensional . . . . . . . . . . . . . . . 34 3.4.2 Resultats num`erics pel model unidimensional . . . . . . . . . . . . . . . . 35 4 1r m`etode de control 38 4.1 Cas d’una sola c`el·lula ................................. 38 4.1.1 Estabilitat del m`etode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 7 Figura 1.4: Propagaci´o normal del potencial d’acci´o en teixit (a), ona espiral associada a la taquicardia (b), trencament de les ones espirals associat a la fibril·laci´o (c). Imatge extreta de l’article [1]. del calci intracel·lular. Experimentalment, observem que aquests fen`omens es produeixen quan el per´ıode (T) entre batec i batec ´es curt, ´es a dir, quan augmenta la nostra pulsaci´o. ´ Es per aix`o que, per confirmar l’exist`encia dels alternans, es realitza una prova no invasiva d’esfor¸c, mesurant l’amplitud de l’ona T per a diferents cicles quan s’est`a entre 90 i 110 batecs per minut i detectant si hi ha difer`encies, com les que podem observar a la Fig.1.5. A la llarga, aquestes alteracions poden ser molt greus, ja que produeixen el trencament de les ones espirals que es donen en el cas de les taquic`ardies, suposant aix´ı el pas a la fibril·laci´o i, si no s’atura el cor amb una desc`arrega el`ectrica mitjan¸cant un desfibril·lador, a la mort. En aquest treball, ens centrarem en els alternans produ¨ıts per variacions en el potencial d’acci´o. ´ Es per aquest motiu que introduirem els termes APD o Action Potencial Duration/Durada del potencial d’acci´o i DI o Diastolic Interval/Interval diast`olic. En acabar el potencial d’acci´o, existeix un interval de temps previ a la seg¨uent activaci´o el qual anomenem di`astole. En aquest interval els ventricles s’omplen de sang despr´es de la contracci´o. Si aquest ´es massa curt, les c`el·lules no hauran pogut recuperar les propietats el`ectriques corresponents, donant a lloc a una durada del potencial d’acci´o m´es curta. Vindr`a seguit d’una di`astole m´es llarga que seguir`a un nou potencial d’acci´o m´es llarg, suposant aix´ı una variaci´o del cicle del potencial d’acci´o. 14 Figura 1.5: Alternans a l’electrocardiograma. Imatge extreta de l’article [2] Actualment, els metges obten per la implantaci´o de desfibril·ladors en pacients amb alternans, ja que com hem exposat abans, les seves conseq¨u`encies a la llarga s´on molt greus. Tamb´e existeixen diferents m`etodes de control per a corregir aquest fenomen, introdu¨ıbles en un dispositiu del tipus marcapassos. ´ Es en l’estudi de diferents m`etodes en el que se centra aquest treball. 1.4 Objectius del treball Vista la base bi`ologica en la que es basa el treball, centrarem l’estudi i posterior correcci´o dels alternans en l’acompliment dels seg¨uents objectius: •A partir d’un model simplificat del comportament el`ectric d’una c`el·lula card´ıaca, estudiar el sistema din`amic resultant i les seves inestabilitats (les quals seran els alternans). •La resoluci´o dels models del comportament el`ectric d’una sola c`el·lula i del teixit unidimensional resultant de considerar el seu acoblament, mitjan¸cant m`etodes num`erics per equacions diferencials ordin`aries i per equacions amb derivades parcials. •L’aplicaci´o de dos m`etodes de control discrets del per´ıode Tper tal de corregir els alternans que apareixen en els dos models. •La comparaci´o i l’extracci´o de conclusions respecte l’efici`encia d’aquest m`etodes en tots dos models. 15 Cap´ıtol 2 Model matem`atic En aquest treball, tractarem els alternans tant des de la perspectiva d’una sola c`el·lula, com des d’un teixit unidimensional d’una certa longitud L, el qual descriur`a el comportament que segueix un seguit de c`el·lules tenint en compte els seus enlla¸cos. 2.1 Model d’una sola c`el·lula (o zero dimensional) La membrana cel·lular ´es semipermeable, permetent el transport d’ions a trav´es de canals i`onics. Aix`o suposa que hi hagi una difer`encia de concentracions i`oniques dins i fora la c`el·lula, i en conseq¨u`encia, una difer`encia de potencial. Podem interpretar que la membrana cel·lular es comporta com un condensador que emmagatzema c`arrega dels diferents ions que la atravessen. ´ Es a dir, segueix el model descrit a l’esquema el`ectric de la Fig 2.1. Les dos corrents que hem de tenir en compte s´on: la Iion, que ´es deguda a les c`arregues que travessen la membrana cel·lular; i la Ic, deguda a l’emmagatzematge de c`arrega a la membrana. Donat que la c`arrega total es conserva, la seva suma ha de ser zero: Ic+Iion = 0 (2.1) Per altra banda la c`arrega que emmagatzema la membrana ´es Q=CmVm, on Cm´es la capacit`ancia i el corrent ´es el canvi en la c`arrega per unitat de temps, per tant Ic=dQ dt =Cm dVm dt (2.2) 16 Figura 2.1: Esquema el`ectric del potencial transmembrana. Imatge extreta de l’article [1] Per tant, l’equaci´o del circuit ´es: Cm dVm dt +Iion = 0 (2.3) amb Vm=Vi−Veque descriuen els voltatges a l’interior i l’exterior de la c`el·lula respectivament. En general podrem descriure els corrents i`onics utilitzant l’expresi´o IX=gX(Vm−VNernst,X ) on Xser`a l’i´o que considerem (Na2+,Ca2+,K+,...). VXcorrespondr`a al potencial de Nernst de l’i´o XigX= 1/rXser`a la conduct`ancia de la membrana per aquest mateix. Com s’exposa a l’ap`endix ??, el potencial de Nernst ´es equivalent a: VNernst,X =Vi,X −Ve,X =kBT qln(ce,X/ci,X) (2.4) on Vi,X iVe,X es corresponen amb els voltatges interior i exterior def ionguts a les diferents concentracions dels ions del tipus Xrespectivament, ci,X ice,X s´on les concentracions d’aquest al interior i exterior de la c`el·lula, kBcorrespon a la constant de Bolzmann, Ta la temperatura iq=z|e|a la c`arrega de l’i´o (z´es la val`encia i |e|la c`arrega d’un electr´o). El primer model d’aquest tipus que es va plantejar va ser el model de Hodgkin-Huxley l’any 1952, el qual descrivia el comportament de les neurones, on el corrent i`onic era: Iion =gNa(Vm−VNa) + gK(Vm−VK) + gL(Vm−VL) (2.5) on IL´es l’anomenada l.leak corrent”. Posteriorment, seguint amb la mateixa formulaci´o, es van desenvolupar altres models com el de Noble [6] l’any 1962 i el de Beeler-Reuter [5] al 1977. Actualment, existeixen molts models que inclouen diferents tipus de c`el·lules card´ıaques que inclouen les humanes. En podem trobar diferents exemples a [8]. 17 Figura 2.2: Potencial d’acci´o de rata. Imatge extreta de l’article [3] En aquest treball considerarem un model simplificat de comportament dels mi`ocits en el que nom´es hi actuen ions de potassi i sodi el qual s’acostuma a donar en animals petits com les rates. Aquesta mancan¸ca de calci d´ona a lloc a qu`e els potencials d’acci´o siguin molt m´es triangulars com els que podem observar a la figura Fig. 2.2. Aquest model l’extraiem de l’article [9]. De manera que el potencial d’acci´o vindr`a descrit per l’equaci´o diferencial: dV dt =−INa +IK+Istim Cm (2.6) On les corrents INa iIKes corresponen als corrents del sodi i el potassi respectivament i Istim ser`a l’est´ımul que introduirem per tal d’excitar el sistema, el qual correspon a l’impuls donat pel node sinusal. Els corrents que pendrem dependran de S(V), que es correspon a una funci´o de Heaviside suavitzada com podem observar a la Fig. 2.3, S(V) = 1 + tanh V−Vc  2−−→ →0S(V) =    1V > VC 0V < VC (2.7) El corrent del potassi vindr`a descrit per l’Eq. (2.8), la qual correspondr`a a una recta creixent per V < VCi una funci´o constant per V > VC, com podem observar a la Fig. 2.4: f IK=IK Cm =1 τ0S(V) + (1 −S(V)) V VC−−→ →0f IK(V) =    1 τ0V > VC 1 τ0 V VCV < VC (2.8) 18 Figura 2.3: Funci´o de Heaviside suavitzada que utilitzem S(V) Figura 2.4: Representacions gr`afiques de les corrents de Potassi i Sodi per h= 1 En el cas del corrent del sodi, INa, actua en sentit contrari i dep`en de la probabilitat que el canal de sodi estigui obert, donada per la variable h(V). La seva equaci´o ser`a doncs l’Eq. (2.9), i suposant h= 1, ser`a de nou una funci´o esglaonada suavitzada com podem veure a Fig. 2.4. g INa=INa Cm =−S(V)h τA −−→ →0g INa(V) =    −h/τAV > VC 0V < VC (2.9) Pel que fa a la probabilitat h, tamb´e ve descrita per una equaci´o ordin`aria en termes de Vi S(V). dh dt =1−S(V)−h τ−(1 −S(V)) + τ+S(V)−−→ →0 dh dt (V) =    −h/τ+V > VC (1 −h)/τ−V < VC (2.10) Per ´ultim considerarem un corrent d’estimulaci´o del sistema Istim que ´unicament dependr`a del temps. Aquest corrent vindr`a descrit per una funci´o pols T-peri`odica com podem observar a la 19 Figura 2.5: Corrent d’estimulaci´o en funci´o del temps par`ametre VCτ0τAτ+τ−H tt valor 0.1 150 6 12 60 0.005 0.015 10 Taula 2.1: Taula de valors dels par`ametres (Fig. 2.5). ] Istim(t) = Istim Cm (t) =    Ht (mod T) < tt 0 altrament (2.11) on H´es la intensitat d’aquesta estimulaci´o, tt ´es l’amplitud de l’interval de temps en qu`e apliquem l’estimulaci´o i T´es el per´ıode que defineix cada quant l’apliquem. De manera que com a resultat tenim el sistema d’equacions ordin`aries seg¨uent: dV dt =−IK+INa +Istim Cm =−S(V) + [1 −S(V)]V/Vc τ0 +S(V)h τA +Istim(t) (2.12) dh dt =1−S(V)−h τ−[1 −S(V)] + Sτ+ (2.13) El qual, per →0, t´e la forma: dV dt =   h/τa−1/τ0V > VC −V/VCτ0V < VC (2.14) dh dt =   −h/τ+V > VC (1 −h)/τ−V < VC (2.15) Per altra banda, prendrem els valors descrits a la taula 2.1, per la implementaci´o dels m`etodes, on les unitats de les variables de temps s´on mil·lisegons. 20 Figura 2.6: Prova num`erica de l’excitabilitat del sistema, resultats obtinguts imposant les condicions inicials V0amb valors entre 0.2 i 1.4 i h0= 1. 2.1.1 Sistema excitable Com veiem a [16], a la biologia ´es prou corrent l’aparici´o de sistemes excitables. Aquests es defineixen per tindre un punt fix, al que si se li apliquen petites pertorvacions torna r`apidament, per`o en el cas que aquestes superin un valor llindar la resposta i la tornada al punt fix ´es molt m´es extensa. Les c`el·lules card´ıaques en s´on un clar exemple, en les que quan es supera VC, degut a l’entrada d’ions de sodi es genera un fort augment del potencial que posteriorment decau amb la sortida dels ions de potassi, a aquesta resposta l’anomenem potencial d’acci´o. Com podem veure a la Fig. 2.6, per valors menors a VCtornen r`apidament a 0 que ´es el punt fix del sistema, en canvi quan arribem a aquest potencial llindar, la resposta del sistema ´es qualitativament molt diferent, trigant un temps en tornar al punt d’equilibri. 2.1.2 Resoluci´o num`erica Per a resoldre l’equaci´o diferencial plantejada a l’Eq. (2.13), discretitzarem respecte al temps considerant un cert ∆tconstant i usarem el m`etode de Dormand-Prince, un m`etode expl´ıcit de 21 la fam´ılia dels Runge-Kutta d’ordre 5. La seva matriu de Butcher ´es: 0 1 5 1 5 3 10 3 40 9 40 4 5 44 45 −56 15 32 9 8 9 19372 6561 −25360 2187 64448 6561 −212 729 19017 3168 −355 33 46732 5247 49 176 −5103 18656 135 384 0500 1113 125 192 −2187 6784 11 84 De manera que amb la notaci´o ˙z(t) = ( ˙ V(t),˙ h(t)) = g(t, (V(t), h(t))), zi= (V(ti), h(ti)) = (V(i·∆t), h(i·∆t)), calculem la soluci´o aplicant l’algorisme seg¨uent: Partint de z0= (V0, h0), condicions inicials donades, Per i= 0,1,2, ..., npassos k1=g(ti, zi) k2=g(ti+1 5∆t, zi+ ∆t·(1 5k1)) k3=g(ti+3 10∆t, zi+ ∆t·(3 40k1+9 40k2)) k4=g(ti+4 5∆t, zi+ ∆t·(44 45k1−56 15k2+32 9k3)) k5=g(ti+8 9∆t, zi+ ∆t·(19372 6561 k1−25360 2187 k2+64448 6561 k3−212 729k4)) k6=g(ti+ ∆t, zi+ ∆t·(9017 3168k1−355 33 k2+46732 5247 k3+49 176k4−5103 18656k5)) zi+1 =zi+ ∆t·(35 384k1+500 1113k3+125 192k4−2187 6784k5+11 84k6)) Seg¨uent iteraci´o. Aquest m`etode ´es el que implementa el Matlab a la funci´o ode45 juntament amb un modificador de pas. Nosaltres no usem aquesta funci´o directament, ja que prenent el pas constant, ens assegurem passar per instants de temps que comptin amb l’impuls Istim. Prenent diferents valors de T, interval que indica cada quant hi ha est´ımuls per part de Istim, observem que el sistema es comporta de maneres molt diferents. Per a Tprou gran, observem que la funci´o V(t) ´es T-peri`odica, i que les ones tenen totes la mateixa al¸cada. En canvi, si T ´es prou petit, el per´ıode de la funci´o V(t) augmenta a 2Tdonant a lloc, dues ones de diferents al¸cades que s’alternen, com podem veure a la Fig. 2.7. 22 Podem trobar el codi a l’ap`endix D.1.1. Figura 2.7: Representaci´o dels resultats obtinguts mitjan¸cant el m`etode de resoluci´o del problema zero dimensional per a T= 400 ms i T= 330 ms. 2.2 Model del teixit card´ıac (unidimensional) En aquest cas, ja no tenim una sola c`el·lula, sin´o un teixit format per un conjunt de c`el·lules que es comporten com les descrites a l’apartat anterior, per tant, tant el voltatge com la probabilitat que els canals de sodi estiguin oberts no dependran nom´es de l’instant de temps, sin´o que tamb´e dependran de la posici´o de cada c`el·lula (V(x, t) i h(x, t)). Donat que coneixem el comportament de cada c`el·lula nom´es cal estudiar el comportament dels acoblaments entre elles. El corrent el`ectric flueix d’una c`el·lula a l’altra a trav´es de les unions de gap. Aquest fenomen, d´ona lloc a la difusi´o del potencial de membrana, que es pot descriure matem`aticament com: ∂V ∂t =∇ · (D∇V)−IIon Cm (2.16) on D´es el coeficient de difusi´o. La qual anomenem equaci´o del cable que s’explica amb m´es detall a l’ap`endix B. De manera que partint de les equacions plantejades a (2.13) afegint el terme de difusi´o ∇ · (D∇V) = D∂2V ∂x2(ja que ens trobem en teixit unidimensional) obtenim: ∂V ∂t =D∂2V ∂x2−S(V) + (1−S(V))V Vc τ0 +S(V)h τA +pols(t, x) (2.17) ∂h ∂t =1−S(V)−h τ−(1 −S(V)) + Sτ+ (2.18) Amb condicions de contorn de Neumann en el voltatge tant a x= 0 i x=L, ´es a dir, ∂V ∂x (0, t) = ∂V ∂x (L, t) = 0. 23 Figura 3.2: Explicaci´o gr`afica del m`etode d’obtenci´o de punts de la corba de restituci´o ´unicament de l’interval diast`olic anterior, com s’explica a l’article [17]. 3.2.2 C`alcul anal´ıtic de la corba de restituci´o Recordem les equacions diferencials que descriuen el comportament el`ectric d’un miocit, s’aproximen per →0 a: Per V > VC:   dV dt =h/τa−1/τ0 dh dt =−h/τ+ (3.3) Per V < VC:   dV dt =−V/VCτ0 dh dt = (1 −h)/τ− (3.4) Considerant la condici´o inicial com a l’inici del potencial d’acci´o (V(0), h(0)) = (VC, h0), tenim que fins a arribar de nou al voltatge cr´ıtic V(t∗) = VC, de fet per t∗=APD, les equacions per les quals es regir`a el sistema s´on les Eq. (3.3). Per altra banda, tamb´e coneixem que des de l’instant de temps t=−DI fins a l’inici, les equacions que descriuran el moviment s´on les (Eq. 3.4). Considerarem les condicions de vora exposades a la taula 3.1, en les que veiem que la probabilitat h(−DI) = h(APD)≈0 ja que durant el potencial d’acci´o ´es exponencialment decreixent, suposant τ+<< AP D. De fet, si no ´es compl´ıs aquesta premissa no podr´ıem afirmar que 30 Figura 3.3: Punts de la corba de restituci´o obtinguts amb l’estabilitzaci´o del sistema, i posterior imposici´o d’un per´ıode. Figura 3.4: Condicions de contorn per tal de trobar l’expressi´o anal´ıtica de la corba de restituci´o. l’APD dep`en ´unicament del DI anterior, sin´o que tamb´e ho faria del potencial d’acci´o anterior, donant a lloc a un sistema d’ordre superior. De manera que l’equaci´o que descriur`a la probabilitat que els canals de sodi estiguin oberts entre els intants t= 0 i t=APD si imposem h(0) = h0ser`a: h(t) = h0e−t/τ+(3.5) Aix´ı que, utilitzant aquesta expressi´o a l’equaci´o diferencial del voltatge tamb´e la podem integrar: dV dt =h0 τa e−t/τ+−1 τ0 ⇒V(t) = ¯ V−τ+h0 τa e−t/τ+−t τ0 (3.6) Tenint en compte V(0) = VC=¯ V−τ+h0/τa⇒¯ V=VC−τ+h0/τa. Per tant, V(t) = VC−τ+h0 τa (1 −e−t/τ+)−t τ0 (3.7) 31 t = -DI t = 0 t = APD VVCVCVC h 0 h00 Taula 3.1: Condicions de vora Per altra banda, tamb´e coneixem V(APD) = VC. Per tant, V(APD) = VC−τ+h0 τa (1 −e−APD/τ+)−APD τ0 =VC ⇒APD =τ0τ+h0 τa (1 −e−APD/τ+) (3.8) Suposant AP D >> τ+, tindrem que e −AP D τ+≈0, aix´ı que APD ≈τ0τ+h0 τa (3.9) Imposant continu¨ıtat, h(0) = h0, i que hexpressa una probabilitat i per tant h≤1, podem integrar-la entre els instants t=−DI it= 0. h(t) = 1 −¯ he−t/τ−⇒h(0) = 1 −¯ h=h0⇒¯ h= 1 −h0(3.10) h(t)=1−(1 −h0)e −t τ−(3.11) Per altra banda tamb´e coneixem h(−DI) = 0, aix´ı que podem trobar h0en funci´o de DI h(−DI)=1−(1 −h0)eDI/τ−= 0 ⇒h0= 1 −e−DI/τ−(3.12) De manera que substitu¨ınt a l’ Eq. (3.9), trobem finalment la relaci´o entre l’APD i el DI anterior. APD ≈τ0τ+ τa (1 −e−DI/τ−) =: f(DI) (3.13) 3.3 Estabilitat del sistema Plantegem doncs l’estudi de l’estabilitat del sistema enunciat a l’apartat de corba de restituci´o, on suposarem que T=DIn+APDn´es fixe: APDn+1 =f(DIn) = f(T−APDn) Per tal d’alleugerir una mica la notaci´o definim an:= AP Dnidn:= DIn. De manera que an+1 =f(dn) = f(T−an) (3.14) 32 Suposem que existeix un punt fix del sistema el qual anomenem a∗=f(T−a∗). Definim la successi´o real zn=an−a∗, de manera que an=a∗+zn. Estudiarem l’estabilitat del sistema descrit a l’Eq. (3.14) en funci´o d’aquesta successi´o zn. an+1 =f(T−an) a∗+zn+1 =f(T−a∗−zn) (3.15) ≈f(T−a∗)−df dd(T−a∗)zn(3.16) ⇒zn+1 =−f0zn(3.17) On de l’Eq. (3.15) a l’Eq. (3.16) hem fet l’expansi´o de Taylor de primer grau al voltant del punt T−a∗i definim f0:= df dd(T−a∗). Per altra banda, experimentalment observem que el fenomen que es produeix, quan apareixen aquestes inestabilitats, ´es que dupliquem el per´ıode de V(t), de manera que an+2 =ani, en conseq¨u`encia, zn+2 =zn∀n∈ N . De manera que, prenent la variable auxiliar yn=zn+1, podem expressar aquest nou sistema com:  yn+1 zn+1  = 0 1 0−f0 = yn zn (3.18) Per tal de fer un estudi te`oric d’aquestes inestabilitats, calcularem el valor del m`odul dels valor propis en funci´o de f0. −λ1 0−f0−λ=λ(f0+λ) = 0 (3.19) Els valors propis seran doncs λ1= 0 i λ2=−f0. Per tant el sistema ser`a estable i convergir`a an−−−→ n→∞ a∗si |λ2|=| − f0|<1 i ser`a inestable i apareixeran alternans si |λ2|=| − f0|>1. Donat que la corba de restituci´o ´es creixent i per tant f0>0, podrem resumir l’estudi de la estabilitat del sistema de l’Eq. (3.14) com: f0<1⇒ESTABLE f0>1⇒INESTABLE Podem observar mitjan¸cant la Fig. 3.5, que iterant per valors prou significatius del per´ıode T, per tot DIinicial, es tendeix al punt fix o equivalentment la intersecci´o de les corbes, com podem veure al primer gr`afic. Per altra banda, si el per´ıode ´es petit, l’efecte ´es el contrari, allunyant-nos del punt fix. 33 Figura 3.5: Convergencia del m`etode en funci´o de f0 3.4 Resultats num`erics Per tal de realitzar un estudi num`eric dels efectes dels alternans, iterarem el programa de resoluci´o del model zero dimensional per a diferents per´ıodes, de fet per a per´ıodes entre 320 i 380 ms, per a veure l’efecte que t´e en les variables: APD, DI i ∆Pics. 3.4.1 Resultats num`erics pel model zero dimensional Iterant el programa de resoluci´o del model zero dimensional per a per´ıodes entre 320 i 380 ms, obtenim els resultats de la Fig. 3.6. En aquests podem observar que no apareixen alternans per a per´ıodes majors a 344ms, on apareix una bifurcaci´o de Pitchfork. De manera que, tot i existir una soluci´o per la que no apar`eixen alternans, aquesta ´es inestable i per aix`o no apareix en les simulacions. ´ Es a partir del m`etodes de control que estabilitzarem aquesta soluci´o suposant aix´ı la correcci´o dels alternans. Per tal de comparar aquest resultat amb la corba de restituci´o aproximada a l’ap`endix C, calculem el valor pel qual f0= 1: APD =f(DI) = a(eb·DI −1) ⇒f0=abeb·DI (3.20) f0=abeb·DI = 1 ⇒DI =1 bln(1/ab) = 94.3348 (3.21) 34 Figura 3.6: Estudi dels efectes dels alternans en la difer`encia entre dos pics consecutius, la durada del potencial d’acci´o i de l’interval diast`olic en funci´o del per´ıode inicial Tpel problema zero dimensional. Per altra banda, substitu¨ınt aquest valor a la corba de restituci´o, tenim: APD =a(eb·DI −1) = 244.3176 ⇒T=APD +DI = 338.6524 (3.22) Aquest valor ´es prou proper a T= 344 ms, fet que ens permet confirmar que aquesta corba est`a ben definida. 3.4.2 Resultats num`erics pel model unidimensional Per la implementaci´o num`erica del model del teixit no hi ha prou amb tenir en compte les dues ´ultimes iteracions, ja que els alternans apareixen en forma de modulacions. Per tant, definirem d’una altra manera la variable ∆Pics. 35 280 290 300 310 320 0 0.2 0.4 Pics(T) al principi per a L = 10 280 290 300 310 320 0 0.2 0.4 Pics(T) al mig per a L = 10 280 290 300 310 320 0 0.2 0.4 Pics(T) al final per a L = 10 0 5 10 15 20 0 0.5 1 Pics(L) al principi per a T = 285 0 5 10 15 20 0 0.5 1 Pics(L) al mig per a T = 285 0 5 10 15 20 0 0.5 1 Pics(L) al final per a T = 285 Figura 3.7: Estudi de les inestabilitats en funci´o del per´ıode inicial per L= 10 cm i en funci´o de la longitud del teixit per el per´ıode inicial T= 285 ms. En primer lloc, calcularem el valor mitj`a de les al¸cades dels ´ultims pics. El nombre de pics variar`a en funci´o del nombre d’iteracions de control o d’estabilitzaci´o que plantejem (al codi considerem l’´ultim quart d’aquestes). Un cop conegut aquest valor mitj`a, calcularem el promig de les difer`encies entre els Pics i aquesta mitjana. De manera que ∆Pics queda definit com: Picmitja =PnP ics i=1 Pici nPics (3.23) ∆Pics := PnP ics i=1 |Pici−Picmitja| nPics =PnP ics i=1 |Pici−PnP ics i=1 Pici nP ics | nPics (3.24) En el cas del model unidimensional, hi haur`a dos factors a tenir en compte alhora d’estudiar els alternans; la longitud del teixit Li el per´ıode inicial T. De manera que cada cop que plantegem l’estudi d’un m`etode de control, fixarem una d’aquestes variables. Pel que fa a la variable espacial, tamb´e ser`a un factor en l’amplitud dels alternans el punt del teixit on ens trobem, aix´ı que sempre els mesurarem a l’inici del teixit, al punt mig d’aquest i al punt final, com podem observar a la Fig. 3.7. En aquest cas, el model ´es estable per valors de Tmajors a 320 ms. Aquesta cota ´es molt menor a la del model unicel·lular, on era estable per valors majors a T= 344 m,s per tant haurem de 36 centrar l’estudi per valors menors als que plantejavem al cas zero dimensional (on anavem de T= 320 ms a T= 380 ms), com per exemple de T= 280 ms a T= 320 ms. 37 Cap´ıtol 4 1r m`etode de control Un cop estabilitzats els alternans per un Tinici fix, controlarem la difer`encia de voltatge usant el per´ıode com a variable discreta Tn, la qual definirem com: Tn=τ+γ 2(APDn−APDn−1) (4.1) on γ´es una constant, τ=Tinici, i AP DniAPDn−1es corresponen als APDs de les iteracions n-`essima i (n-1)-`essima respectivament. A simple vista, aquest m`etode de control augmenta el temps entre pulsacions en cas que l’anterior APD fos m´es curt que l’actual i el disminueix en cas contrari. Podem trobar aplicacions d’aquest m`etode de control a [10] i [11]. 4.1 Cas d’una sola c`el·lula En aquesta secci´o ens centrarem en el primer model. Per come¸car, delimitarem els valors de γpels quals el m`etode de control funciona. Un cop fixats uns criteris sobre aquesta variable, observarem quins resultats num`erics t´e al aplicar-los: tant per a un per´ıode fix, com per delimitar el rang de per´ıodes pel qual el control funciona. Podem trobar el codi a l’ap`endix D.2.1. 38 4.1.1 Estabilitat del m`etode De nou, considerarem que existeix a∗punt fix del mapa an+1 =f(dn) = f(Tn−an), ´es a dir, a∗=f(d∗) = f(τ−a∗). Definirem ancom l’AP D a la n-`essima iteraci´o, el qual podrem descriure com an=a∗+znon {zn}n⩾0´es una successi´o discreta. Prenent el Tndefinit a (4.1) tenim: an+2 =a∗+zn+2 =f(Tn+1 −an+1) (4.2) =f([τ+γ 2([a∗+zn+1]−[a∗+zn])] −[a∗+zn+1]) (4.3) =f(τ−a∗) + f0[γ 2(zn+1 −zn)−zn+1] (4.4) ⇒zn+2 =f0[γ 2(zn+1 −zn)−zn+1] (4.5) On f0:= ∂f ∂d (T−a∗). Prenent la variable auxiliar yn=zn+1 tenim el sistema:    yn+1 =f0[γ 2(yn−zn)−yn] zn+1 =yn=⇒ yn+1 zn+1  = f0(γ 2−1) f0γ 2 1 0   yn zn  Que tindr`a per valors pr`opis, les solucions de l’equaci´o:  f0(γ 2−1) −λ f0γ 2 1−λ=λ2−f0(γ 2−1)λ+f0γ 2= 0 Recordem que el m`etode ser`a estable si |λ|<1 i inestable si |λ|>1. Resolent la equaci´o de segon grau tenim: λ=1 2[f0(γ 2−1) ±r(f0)2(γ 2−1)2−4f0γ 2] (4.6) Donat que per a valors de f0<1 a l’iterar amb el mateix per´ıode inicial no apareixen inestabilitats i en conseq¨u`encia no hi haur`a alternans, ens centrarem en el cas f0>1. El nostre objectiu ser`a acotar γ(f0) per tal de que el control sigui efectiu. Partint de γ= 0, observem que els vaps s´on λ= 0 i λ=−f0, per tant busquem el m´ınim γtal que tots dos VAPS siguin menors en m`odul a 1. •Si q(f0)2(γ 2−1)2−4f0γ 2∈R⇒λ∈R, per l’exposat anteriorment ens centrarem en el valor frontera d’estabilitat, ´es a dir, λ=−1. Si λ=−1, (−1)2−f0(γ 2−1)(−1) + f0γ 2= 0 ⇒1 + f0(γ−1) = 0 ⇒γ= 1 −1 f0. Per tant, la cota inferior de γper la qual el m`etode ser`a estable ´es: γ > 1−1 f0(4.7) 39 280 290 300 310 320 0 0.2 0.4 Pics(T) al principi per a L = 10 pre Control APD Control 280 290 300 310 320 0 0.2 0.4 Pics(T) al mig per a L = 10 280 290 300 310 320 0 0.5 Pics(T) al final per a L = 10 0 5 10 15 20 0 0.5 1 Pics(L) al principi per a T = 285 pre Control APD Control 0 5 10 15 20 0 0.5 1 Pics(L) al mig per a T = 285 0 5 10 15 20 0 0.5 1 Pics(L) al final per a T = 285 Figura 4.6: Resultats de l’estudi de les influ`encies del per´ıode inicial i la longitud del teixit en el control de l’APD aplicat al model unidimensional. 46 Cap´ıtol 5 2n m`etode de control, DI constant Novament, centrarem el control dels alternans prenent el per´ıode entre els pols el`ectrics com una variable discreta, la qual modificarem a cada iteraci´o. Definirem doncs el per´ıode a la n-`essima iteraci´o com: Tn=APDn+DI∗(5.1) On el valor DI∗=: d∗´es un valor fixat. De manera que aquest m`etode consistir`a en la imposici´o d’un interval diast`olic concret, amb el qual esperem obtenir una resposta constant en la durada del potencial d’acci´o, ja que coneixem que nom´es dep`en d’aquest DI, com hem vist a l’apartat d’estabilitat. Podem trobar aplicacions d’aquest m`etode de control a [12] i [13]. 5.1 Cas d’una sola c`el·lula En primer lloc ens centrarem en el primer model, ´es a dir, en el d’una sola c`el·lula. Per aquest model, estudiarem l’estabilitat del m`etode per tal de provar la seva efic`acia. En segon lloc, parlarem de com implementar-lo num`ericament, ja que com observarem a la secci´o 5.1.1, per tal que funcioni aquesta variable ha de ser exactament d∗=T−a∗i per tant necessitem coneixements previs de la corba de restituci´o. Finalment, estudiarem els resultats que obtenim amb aquesta implementaci´o. Podem trobar el codi a l’ap`endix D.3.1. 47 5.1.1 Estabilitat aplicant el m`etode Considerem de nou el sistema APDn+1 =f(DIn) o equivalentment, an+1 =f(dn) = f(T−an) que t´e un punt fix a a∗=f(T−a∗). De manera que tenim an+1 =f(Tn−an) (5.2) =f(an+d∗−an) (5.3) =f(d∗) (5.4) Per tant, prenent d∗=T−a∗assegurem la converg`encia del m`etode. 5.1.2 Implementaci´o num`erica Per la implementaci´o num`erica d’aquest m`etode necessitem coneixements previs del sistema per tal de poder imposar el d∗adient, sense alterar el per´ıode de batec. En concret, hem de con`eixer l’expressi´o de la corba de restituci´o. A l’annex C expliquem com a partir de valors coneguts en podem calcular una expressi´o, utilitzant el m`etode dels m´ınims quadrats. Un cop definida la nostra aproximaci´o de la corba de restituci´o, per tal d’imposar el d∗adient per a cada per´ıode T, buscarem la intersecci´o de la recta T=AP D +DI iAPD =f(DI) calculant el zero de g(DI) = f(DI)−(T−DI). Per fer-ho, usarem la funci´o de Matlab fsolve, utilitzant com a punt inicial un promig dels DI’s. Un altre aspecte a tenir en compte en la implementaci´o num`erica d’aquest control ´es que el DI d’una iteraci´o inclou el temps de reacci´o un cop activat l’impuls el`ectric. Per tal de tindre’l en compte, un cop hagem iterat la meitat de per´ıodes de control, restarem aquest temps de reacci´o al T que imposem. 5.1.3 Resultats num`erics per un per´ıode inicial determinat Novament, prenent un per´ıode per al qual hi ha alternans en un inici com T= 330 ms i aplicant el control, les inestabilitats es corregeixen com podem observar a la (Fig. 5.1). En aquesta figura hem aplicat 20 iteracions de control despr´es d’haver estabilitzat el sistema durant 20 iteracions m´es. Contemplem que abans de corregir el DI tenint en compte el temps de reacci´o, el per´ıode tendia a ser una mica m´es alt de l’inicial, per`o un cop aplicada la correcci´o, tendeix a aquest. 48 Figura 5.1: Resultats per T= 330 ms en el cas d’una sola c`el·lula 5.1.4 Resultats num`eric en funci´o del per´ıode inicial Pel que fa a l’efectivitat del m`etode en funci´o del per´ıode inicial, igual que en el primer control, aquest funciona per tot per´ıode inicial que considerem, com podem observar a la Fig. 5.2. 5.2 Cas en teixit De nou en el cas unidimensional, aplicarem els resultats obtinguts del model d’una sola c`el·lula. Com hem fet en el cas del control anterior, tamb´e tindrem en compte els factors de la longitud del teixit i el per´ıode inicial per tal d’estudiar l’efectivitat del m`etode de control. Podem trobar el codi a l’ap`endix D.3.2. 5.2.1 Resultats num`erics per un per´ıode i una longitud determinades Prenent T= 290, per´ıode inicial pel qual sabem que apareixen alternans en el cas unidimensional amb L= 10, i aplicant el control, observem que els alternans es corregeixen a l’inici del teixit, per`o segueixen al mig i al final d’aquest com podem veure a la Fig. 5.3, on hem aplicat 20 iteracions d’estabilitzaci´o i 20 de control. Com en el primer control, si mantenim la longitud del teixit per`o prenem un per´ıode m´es gran, en aquest cas T= 305, el control s’est´en a tot el teixit eliminant els alternans d’arreu com podem observar a la Fig. 5.4 on hem aplicat el mateix nombre d’iteracions. 49 Figura 5.2: Estudi dels resultats del segon m`etode de control pel primer model en funci´o del per´ıode inicial Tamb´e, si mantenim el per´ıode T= 290 ms per`o redu¨ım la longitud del teixit en L= 4 cm, aconseguim corregir les inestabilitats com podem observar a la Fig. ?? obtinguda amb el mateix procediment. 5.2.2 Resultats num`erics en funci´o del per´ıode per una longitud determinada Si fixem la variable longitud en L= 10 cm, el rang de valors pels quals el control funciona ´es prou semblant al del primer m`etode com podem observar a la Fig. 5.6. En aquest cas, la correcci´o dels alternans a l’inici del teixit ´es possible sigui quin sigui el per´ıode inicial. Pel que fa a la correcci´o en la totalitat del teixit per`o, aquesta nom´es es produeix per per´ıodes inicials majors a T= 305 ms. 50 0 2000 4000 6000 8000 10000 t 0 1 2Vinici(t), Tinicial =290 i L =10 0 2000 4000 6000 8000 10000 t 0 1 2Vmig(t), Tinicial =290 i L =10 0 2000 4000 6000 8000 10000 t 0 1 2Vfi(t), Tinicial =290 i L =10 0 5 10 15 20 iteració 275 280 285 290 295 300 305 T T =290 i L =10 Figura 5.3: Resultats per T= 290 ms i L= 10 cm 5.2.3 Resultats num`erics en funci´o de la longitud per a un per´ıode determinat Per tal d’observar l’efecte que t´e la longitud en l’efic`acia d’aquest segon m`etode, prenem el per´ıode inicial T= 285 ms, per´ıode pel qual sabem que es produeixen alternans, i resolem per a longituds entre 0.5 i 20 cm, com pr`eviament hav´ıem fet en el m`etode de control anterior. En aquest cas, obtenim resultats pr`acticament id`entics que en el anterior m`etode com veiem a la Fig. 5.6, ja que el control nom´es es d´ona en la totalitat del teixit per longituds menors o iguals aL= 5 cm, i sempre es corregeixen els alternans a l’inici del teixit. 51 0 5000 10000 t 0 1 2Vinici(t), Tinicial =305 i L =10 0 5000 10000 t 0 1 2Vmig(t), Tinicial =305 i L =10 0 5000 10000 t 0 1 2Vfi(t), Tinicial =305 i L =10 0 5 10 15 20 iteració 295 300 305 310 315 320 T T =305 i L =10 Figura 5.4: Resultats per T= 305 ms i L= 10 cm. 0 2000 4000 6000 8000 10000 t 0 1 2Vinici(t), Tinicial =290 i L =4 0 2000 4000 6000 8000 10000 t 0 1 2Vmig(t), Tinicial =290 i L =4 0 2000 4000 6000 8000 10000 t 0 1 2Vfi(t), Tinicial =290 i L =4 0 5 10 15 20 iteració 280 290 300 310 T T =290 i L =4 Figura 5.5: Resultats per T= 290 i L= 4 cm 52 280 290 300 310 320 0 0.2 0.4 Pics(T) al principi per a L = 10 pre Control Constant DI Control 280 290 300 310 320 0 0.2 0.4 Pics(T) al mig per a L = 10 280 290 300 310 320 0 0.5 Pics(T) al final per a L = 10 0 5 10 15 20 0 0.5 1 Pics(L) al principi per a T = 285 pre Control Constant DI Control 0 5 10 15 20 0 0.5 1 Pics(L) al mig per a T = 285 0 5 10 15 20 0 0.5 1 Pics(L) al final per a T = 285 Figura 5.6: Resultats de l’estudi de les influ`encies del per´ıode inicial i la longitud del teixit en el control del DI constant aplicat al model unidimensional. 53 Cap´ıtol 6 Comparaci´o dels m`etodes Un cop vistos aquests dos m`etodes de control dels alternans, tractarem de comparar-los des del punt de vista dels resultats que hem obtingut, com des dels coneixements previs del sistema que necessitem per tal d’implementar-los. Pel que fa als resultats, podem observar que s´on bastant similars tant pel model d’una sola c`el·lula, com pel model de teixit, on ´es lleugerament millor el m`etode de control del DI constant. Pel que fa al primer, la correcci´o de les inestabilitats es d´ona per tot per´ıode inicial, pel que podem concloure que tots dos m`etodes funcionen. En el cas del model en teixit, aquests controls, corregeixen les inestabilitats a l’inici de la fibra per tot per´ıode inicial, per`o el control es perd quan ens allunyem d’aquest com podem observar a la (Fig. 6.1). Tamb´e, tots dos m`etodes, depenen de la longitud del teixit de manera equivalent. Per tant, a nivell de resultats podem concloure que tots dos m`etodes s´on igual d’efica¸cos. Per altra banda, si ens centrem en els coneixements que requereixen, el m`etode de control del DI constant necessita un coneixement del sistema previ a la seva implementaci´o, ja que no es pot aplicar sense con`eixer l’expressi´o de la corba de restituci´o aproximada, com la del ap`endix C. En canvi, l’aplicaci´o del m`etode de l’APD no requereix cap m´es coneixement que l’interval pel qu`e γel fa estable. 54 280 290 300 310 320 0 0.2 0.4 Pics(T) al principi per a L = 10 pre Control APD Control Constant DI Control 280 290 300 310 320 0 0.2 0.4 Pics(T) al mig per a L = 10 280 290 300 310 320 0 0.5 Pics(T) al final per a L = 10 0 5 10 15 20 0 0.5 1 Pics(L) al principi per a T = 285 pre Control APD Control Constant DI Control 0 5 10 15 20 0 0.5 1 Pics(L) al mig per a T = 285 0 5 10 15 20 0 0.5 1 Pics(L) al final per a T = 285 Figura 6.1: Resultats obtinguts amb el m`etode de control de l’APD amb γ= 0.4 i el m`etode de control del DI constant. 55 Prenent el l´ımit quan dx →0, tenim: Ii(x) = gi ∂Vi ∂x (B.4) Ie(x) = ge ∂Ve ∂x (B.5) It=−∂Ii ∂x =∂Ie ∂x (B.6) De manera que substituint a l’Eq. (B.3) i utilitzant Vi=V+Ve, pCm dV dt +Iion=∂ ∂xgi ∂Vi ∂x =−∂ ∂xge ∂Ve ∂x (B.7) pCm dV dt +Iion=∂ ∂xgi ∂V ∂x +∂ ∂xgi ∂Ve ∂x (B.8) Per altra banda, per l’equaci´o l’Eq. (B.6) sabem que, ∂Ii ∂x +∂Ie ∂x =∂ ∂x(Ii+Ie) = ∂ ∂xgi ∂Vi ∂x +ge ∂Ve ∂x = 0 (B.9) suposant que les conduct`ancies s´on homog`enies. ∂ ∂xgi ∂Vi ∂x +ge ∂Ve ∂x =∂ ∂xgi ∂V ∂x + (gi+ge)∂Ve ∂x = 0 (B.10) que t´e per soluci´o: ge ∂Ve ∂x =−gi ∂Vi ∂x +C, o∂Ve ∂x =−gi gi+ge ∂V ∂x +C gi +ge (B.11) d’on obtenim l’equaci´o unidimensional del cable: pCm ∂V ∂t +Iion=∂ ∂xgige gi+ge ∂V ∂x (B.12) De manera que prenent D=1 pCm gige gi+ge, ∂V ∂t =1 pCm ∂ ∂xgige gi+ge ∂V ∂x −Iion Cm =D∂2V ∂x2−Iion Cm (B.13) ja que suposem que les conduct`ancies no dep´enen de la posici´o. 62 Ap`endix C Aproximaci´o de la corba de restituci´o per m´ınims quadrats Donat que a la secci´o 3.2.2 hem vist que la corba de restituci´o es pot aproximar num`ericament per una funci´o exponencial de la forma f(x) = a·(ebx −1), a partir dels valors (DI, APD) obtinguts a la secci´o 3.2.1, buscarem els par`ametres que minimitzin: ¯g(a, b, c) = ||(f(DIi;a, b)−APDi)||2(C.1) O equivalentment, que minimitzin: g(a, b, c) = n X i=1 (a(·ebDIi−1) −APDi)2(C.2) Per fer-ho buscarem (a∗, b∗) tal que (∂g ∂a ,∂g ∂b )(a∗, b∗) = (0,0) on, ∂g ∂a = 2 n X i=1 ebDIi(a(·ebDIi−1) −APDi) ∂g ∂b = 2a n X i=1 ebDIiDIi(a(·ebDIi−1) −APDi) Per resoldre aquest problema usarem la funci´o fsolve del Matlab prenent com a punt inicial (−τ0τ+ tauA,−1 τ−), ja que ´es la soluci´o anal´ıtica que obtenim suposant APD >>> τ+. Com podem veure a la Fig. C.1, la corba resultant es prou propera als valors coneguts de la corba de restituci´o. 63 Figura C.1: Aproximaci´o per m´ınims quadrats de la corba de restituci´o a partir dels valors (DI, APD) de la secci´o 3.2.1 64 Ap`endix D Codis dels programes Tots els codis estan plantejats per la simulaci´o dels models i l’aplicaci´o dels m`etodes de control en aquest utilitzant el Matlab. D.1 Simulaci´o dels models D.1.1 Model d’una sola c`el·lula function [y, t, APD, DI, APics] = estabilitzar_0D(y0, n_periodes, h, vull_APD, vull_V) global T; tend = n_periodes*T; npassos = tend/h; y = []; y(:, 1) = y0; t =[0: h: tend]; for j = 2:npassos+1 y(:, j) = step_DOPRI45_b1(@f, t(j), y(:, j-1), h); end 65 if (vull_APD) [APD, DI, APics] = APD_DI_APics (y (1, (n_periodes-3)*T/h:end), h); else APD = 0; DI = 0; APics = 0; end if (vull_V) y = y; t = [0: h: tend]; else y=0;t=0; end end D.1.2 Model en teixit function [V, h, t, APD, DI, APics] = estabilitzar_1D(V0, h0, n_periodes, Ax, At, L, vull_APD, vull_V) global T; global Vc; APD = []; DI = []; D = 0.001; V = [V0]; h = [h0]; x = [0:Ax:L]; m = length(x); t = [0:At:n_periodes*T]; tsteps = length(t); for j = 1:(tsteps-1) 66 for i = 2:11 F = f(t(j), [V(i,j), h(i,j)]); % els primers 10 punts incorporen el terme del pols V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end for i = 12:(m+1) F = f_sin(t(j), [V(i,j), h(i,j)]); V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end V(1, j+1) = V(2,j+1); V(m+2, j+1) = V(m+1, j+1); end if (vull_APD) APics = zeros(1,3); Pics_inici = []; k = 1; for i = (3/4)*n_periodes:n_periodes-1 Pics_inici(k) = max(V(1, (T/At)*i+1:(T/At)*(i+1))); k = k+1; end Pic_mitjana_inici = sum(Pics_inici)/(k-1); APics(1) = sum(abs(Pics_inici-Pic_mitjana_inici))/(k-1); DI_inici = []; APD_inici = []; j = ((3/4)*n_periodes - 1)*(T/At) + 1; while (V(1, j) < Vc); j = j+1; end while (V(1, j) >= Vc); j = j+1; end % fi APD for i =1:((1/4)*n_periodes - 1) if (j + 2*(T/At) > length(V(1, :)));break; end 67 k = j; %inici DI while (V(1, j) < Vc && j+1 < length(V(1,:))); j = j+1; end DI_inici(i) = (j-k)*At; k = j; %inici APD while (V(1, j) >= Vc && j+1 < length(V(1,:))); j = j+1; end APD_inici(i) = (j-k)*At; end APD(1,:) = APD_inici; DI (1,:) = DI_inici; tol = 1e-3; jj = 1; while(abs(V(ceil(m/2)+1, jj)) < tol); jj = jj+1; end Pics_mig = []; k = 1; i = (3/4)*n_periodes; while ((T/At)*(i+1) + jj < length(V(ceil(m/2) + 1,:))) Pics_mig(k) = max(V(ceil(m/2)+1, (T/At)*i+jj:(T/At)*(i+1) + jj)); k = k+1; i = i+1; end Pic_mitjana_mig = sum(Pics_mig)/(k-1); APics(2) = sum(abs(Pics_mig-Pic_mitjana_mig))/(k-1); DI_mig = []; APD_mig = []; j = ((3/4)*n_periodes - 1)*(T/At) + jj; while (V(ceil(m/2)+1, j) < Vc ) j = j+1 end while (V(ceil(m/2)+1, j) >= Vc) j = j+1 end % fi APD 68 for i =1:((1/4)*n_periodes - 1) if (j + 2*T/At > length(V(ceil(m/2)+1, :))); break; end k = j; %inici DI while (V(ceil(m/2) + 1, j) < Vc && j+1 < length(V(ceil(m/2) +1 ,:))) j = j+1 end DI_mig(i) = (j-k)*At; k = j; %inici APD while (V(ceil(m/2) + 1, j) >= Vc && j+1 < length(V(ceil(m/2) +1 ,:))) j = j+1 end APD_mig(i) = (j-k)*At; end APD(2,1:length(APD_mig)) = APD_mig; DI (2,1:length(DI_mig)) = DI_mig; jj = 1; while(abs(V(m+1, jj)) < tol); jj = jj+1; end Pics_final = []; k = 1; i = (3/4)*n_periodes; while ((T/At)*(i+1) + jj < length(V(m+1,:))) Pics_final(k) = max(V(m+1, (T/At)*i+jj:(T/At)*(i+1) + jj)); k = k+1; i = i+1; end Pic_mitjana_final = sum(Pics_final)/(k-1); APics(3) = sum(abs(Pics_final-Pic_mitjana_final))/(k-1); DI_final = []; APD_final = []; j = ((3/4)*n_periodes - 1)*(T/At) + jj; while (V(m+1, j) < Vc); j = j+1; end while (V(m+1, j) >= Vc); j = j+1; end % fi APD 69 for i =1:((1/4)*n_periodes - 1) if (j + 2*T/At > length(V(m+1, :))); break; end k = j; %inici DI while (V(m + 1, j) < Vc && j+1 < length(V(m+1, :))); j = j+1; end DI_final(i) = (j-k)*At; k = j; %inici APD while (V(m + 1, j) >= Vc && j+1 < length(V(m+1, :))); j = j+1; end APD_final(i) = (j-k)*At; end APD(3,1:length(APD_final)) = APD_final; DI (3,1:length(DI_final)) = DI_final; else APD = 0; DI = 0; APics = 0; end if (not(vull_V)) V=0;h=0;t=0; end end D.2 Aplicaci´o del Control de l’APD D.2.1 Model d’una sola c`el·lula function [y, t, TT, APD, DI, APics] = ControlAPD_0D_v2(y0, n_periodes, h, gamma, APD_0, vull_APD, vull_V) global T; global tt; tau = T; 70 TT = []; t = []; t(1) = 0; APD_ant = APD_0; y = []; y(:, 1) = y0; Vc = 0.1; ant = 1; for j = 1:n_periodes for k = (ant+1):(ant+floor(tt/h)) y(:, k) = step_DOPRI45_b1(@f_con, h*k, y(:, k-1), h); end i = ant; while (y(1, i) < Vc) i=i+1; end inici_APD = i+1; k = ant+floor(tt/h); while (y(1, k) > Vc) k = k+1; y(:, k) = step_DOPRI45_b1(@f_sin, h*k, y(:, k-1), h); end APD_nou = (k-inici_APD)*h; TT(j) = floor(tau + (gamma/2)*(APD_nou - APD_ant)); for i = (k+1):(ant + ceil(TT(j)/h)) 71 else APD = 0; DI = 0; APics = 0; end if (vull_V) t = [At: At: ant*At]; V = V(:,2:end); h = h(:,2:end); else V=0;h=0;t=0; end end D.3 Aplicaci´o del Control del DI constant D.3.1 Model d’una sola c`el·lula function [y, t, TT, APD, DI, APics] = ControlConstantDI_0D(y0, n_periodes, h, DI, vull_APD, vull_V) global tt; TT = []; ant = 1; y = []; y(:, 1) = y0; for j = 1:ceil(n_periodes/2) for i = (ant+1):(ant+floor(tt/h)) y(:, i) = step_DOPRI45_b1(@f_con, (i-1)*h, y(:, i-1), h); end 78 k = ant+floor(tt/h); while (y(1,k) > 0.1) k = k+1; y(:, k) = step_DOPRI45_b1(@f_sin, (k-1)*h, y(:, k-1), h); end for i = (k+1):(k + (DI/h)) y(:, i) = step_DOPRI45_b1(@f_sin, (i-1)*h, y(:, i-1), h); end TT(j) = (k + floor(DI/h) - ant)*h; ant = k + floor(DI/h); end %medim el temps que triga en arribar a Vc la cellula END = length(y(1,:)); i = END - (TT(end)/h); while (y(1, i) < 0.1) i = i+1; end t_reaccio = (i - (END - (TT(end)/h)))*h; %definim el nou DI DI = DI-t_reaccio; for j = ceil(n_periodes/2)+ 1:n_periodes for i = (ant+1):(ant+floor(tt/h)) y(:, i) = step_DOPRI45_b1(@f_con, (i-1)*h, y(:, i-1), h); end k = ant+floor(tt/h); 79 while (y(1,k) > 0.1) k = k+1; y(:, k) = step_DOPRI45_b1(@f_sin, (k-1)*h, y(:, k-1), h); end for i = (k+1):(k + (DI/h)) y(:, i) = step_DOPRI45_b1(@f_sin, (i-1)*h, y(:, i-1), h); end TT(j) = (k + (DI/h) - ant)*h; ant = k + floor(DI/h); end if (vull_APD) [APD, DI, APics] = APD_DI_APics(y(1, (end-((TT(end) + TT(end-1) + TT(end-2))/h):end)), h); else APD = 0; DI = 0; APics = 0; end if (vull_V) t = [h:h:(ant-1)*h]; y = y(:, 2:end); else y=0;t=0; end end 80 D.3.2 Model en teixit function [V, h, t, TT, APD, DI, APics] = ControlConstantDI_1D_v2(V0, h0, n_periodes, Ax, At, L, DI_est, vull_APD, vull_V) global tt; global Vc; TT = []; V = [V0]; h = [h0]; x = [0:Ax:L]; m = length(x); ant = 1; D = 0.001; APD = []; DI = []; APics = []; for jj = 1:ceil(n_periodes/2) k = ant+floor(tt/At); for j = ant:k for i = 2:11 F = f_con(j*At, [V(i,j), h(i,j)]); % els primers 10 punts incorporen el terme del pols V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end for i = 12:(m+1) F = f_sin(j*At, [V(i,j), h(i,j)]); V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); 81 h(i, j+1) = h(i,j) + At*F(2); end V(1, j+1) = V(2,j+1); h(1, j+1) = h(2, j+1); V(m+2, j+1) = V(m+1, j+1); h(m+2, j+1) = h(m+1, j+1); end j = ant+floor(tt/At); while (V(2,j) > 0.1) for i = 2:(m+1) F = f_sin(j*At, [V(i,j), h(i,j)]); V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end V(1, j+1) = V(2,j+1); h(1, j+1) = h(2, j+1); V(m+2, j+1) = V(m+1, j+1); h(m+2, j+1) = h(m+1, j+1); j=j+1; end for k = j:(j + floor(DI_est/At)) for i = 2:(m+1) F = f_sin(j*At, [V(i,k), h(i,k)]); V(i, k+1) = V(i,k) + D*At/(Ax^2)*(V(i+1,k)-2*V(i,k)+ V(i-1, k))+At*F(1); h(i, k+1) = h(i,k) + At*F(2); end V(1, k+1) = V(2,k+1); h(1, k+1) = h(2, k+1); V(m+2, k+1) = V(m+1, k+1); h(m+2, k+1) = h(m+1, k+1); end TT(jj) = (j + floor(DI_est/At) - ant)*At; ant = j + floor(DI_est/At); 82 end %medim el temps que triga en arribar a Vc la cellula END = length(V(1,:)); i = END - (TT(end)/At); while (V(1, i) < 0.1) i = i+1; end t_reaccio = (i - (END - (TT(end)/At)))*At; %definim el nou DI DI_est = DI_est-t_reaccio; for jj = (ceil(n_periodes/2)+1):n_periodes k = ant+floor(tt/At); for j = ant:k for i = 2:11 F = f_con(j*At, [V(i,j), h(i,j)]); % els primers 10 punts incorporen el terme del pols V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end for i = 12:(m+1) F = f_sin(j*At, [V(i,j), h(i,j)]); V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end V(1, j+1) = V(2,j+1); h(1, j+1) = h(2, j+1); V(m+2, j+1) = V(m+1, j+1); h(m+2, j+1) = h(m+1, j+1); end 83 j = ant+floor(tt/At); while (V(2,j) > 0.1) for i = 2:(m+1) F = f_sin(j*At, [V(i,j), h(i,j)]); V(i, j+1) = V(i,j) + D*At/(Ax^2)*(V(i+1,j)-2*V(i,j)+ V(i-1, j))+At*F(1); h(i, j+1) = h(i,j) + At*F(2); end V(1, j+1) = V(2,j+1); h(1, j+1) = h(2, j+1); V(m+2, j+1) = V(m+1, j+1); h(m+2, j+1) = h(m+1, j+1); j=j+1; end for k = j:(j + floor(DI_est/At)) for i = 2:(m+1) F = f_sin(j*At, [V(i,k), h(i,k)]); V(i, k+1) = V(i,k) + D*At/(Ax^2)*(V(i+1,k)-2*V(i,k)+ V(i-1, k))+At*F(1); h(i, k+1) = h(i,k) + At*F(2); end V(1, k+1) = V(2,k+1); h(1, k+1) = h(2, k+1); V(m+2, k+1) = V(m+1, k+1); h(m+2, k+1) = h(m+1, k+1); end TT(jj) = (j + floor(DI_est/At) - ant)*At; ant = j + floor(DI_est/At); end if (vull_APD) APics = zeros(1,3); 84 Pics_inici = []; ant = ceil(sum(TT(1:(3/4)*n_periodes))/At); for i = 1:((1/4)*n_periodes)-1 Pics_inici(i) = max(V(1, ant+1:ant + TT((3/4)*n_periodes+i)/At)); ant = ant + ceil(TT((3/4)*n_periodes+i)/At); end Pic_mitjana_inici = sum(Pics_inici)/length(Pics_inici); APics(1) = sum(abs(Pics_inici-Pic_mitjana_inici))/length(Pics_inici); DI_inici = []; APD_inici = []; j = ceil(sum(TT(1:(3/4)*n_periodes))/At+1); while (V(1, j) < Vc); j = j+1; end while (V(1, j) >= Vc); j = j+1; end % fi APD for i =1:((1/4)*n_periodes - 1) if (j + ceil((TT(end) + TT(end-1))/At) > length(V(1, :))); break; end k = j; %inici DI while (V(1, j) < Vc && j+1 < length(V(1, :))); j = j+1; end DI_inici(i) = (j-k)*At; k = j; %inici APD while (V(1, j) >= Vc && j+1 < length(V(1, :))); j = j+1; end APD_inici(i) = (j-k)*At; end APD(1,:) = APD_inici; DI (1,:) = DI_inici; tol = 1e-3; jj = 1; while(abs(V(ceil(m/2)+1, jj)) < tol); jj = jj+1; end ant = ceil(sum(TT(1:(3/4)*n_periodes))/At) + jj; Pics_mig = []; k = 1; for i = 1:(1/4)*n_periodes-1 85 if(ant + ceil(TT((3/4)*n_periodes+i)/At) > length(V(ceil(m/2) + 1, :))) break end Pics_mig(i) = max(V(ceil(m/2) + 1, ant+1:ant + ceil(TT((3/4)*n_periodes+i)/At))); ant = ant + ceil(TT((3/4)*n_periodes+i)/At); end Pic_mitjana_mig = sum(Pics_mig)/length(Pics_mig); APics(2) = sum(abs(Pics_mig-Pic_mitjana_mig))/length(Pics_mig); DI_mig = []; APD_mig = []; j =ceil(sum(TT(1:(3/4)*n_periodes))/At) + jj; while (V(ceil(m/2)+1, j) < Vc); j = j+1; end while (V(ceil(m/2)+1, j) >= Vc); j = j+1; end % fi APD for i =1:((1/4)*n_periodes - 1) if (j + (TT(end) + TT(end-1))/At > length(V(ceil(m/2)+1, :))); break; end k = j; %inici DI while (V(ceil(m/2) + 1, j) < Vc && j+1 < length(V(1, :))); j = j+1; end DI_mig(i) = (j-k)*At; k = j; %inici APD while (V(ceil(m/2) + 1, j) >= Vc && j+1 < length(V(1, :))); j = j+1; end APD_mig(i) = (j-k)*At; end APD(2,1:length(APD_mig)) = APD_mig; DI (2,1:length(DI_mig)) = DI_mig; jj = 1; while(abs(V(m+1, jj)) < tol); jj = jj+1; end Pics_final = []; k = 1; ant = ceil(sum(TT(1:(3/4)*n_periodes))/At) + jj; for i = 1:(1/4)*n_periodes-1 86 if(ant + ceil(TT((3/4)*n_periodes+i)/At) > length(V(m + 1, :))); break; end Pics_final(i) = max(V(m+1, ant+1:ant + ceil(TT((3/4)*n_periodes+i)/At))); ant = ant + ceil(TT((3/4)*n_periodes+i)/At); end Pic_mitjana_final = sum(Pics_final)/length(Pics_final); APics(3) = sum(abs(Pics_final-Pic_mitjana_final))/length(Pics_final); DI_final = []; APD_final = []; j = ceil(sum(TT(1:(3/4)*n_periodes))/At) + jj; while (V(m+1, j) < Vc); j = j+1; end while (V(m+1, j) >= Vc); j = j+1; end % fi APD for i =1:((1/4)*n_periodes - 1) if (j + (TT(end) + TT(end-1))/At > length(V(m+1, :))); break; end k = j; %inici DI while (V(m + 1, j) < Vc && j+1 < length(V(1, :))); j = j+1; end DI_final(i) = (j-k)*At; k = j; %inici APD while (V(m + 1, j) >= Vc && j+1 < length(V(1, :))); j = j+1; end APD_final(i) = (j-k)*At; end APD(3,1:length(APD_final)) = APD_final; DI (3,1:length(DI_final)) = DI_final; else APD = 0; DI = 0; APics = 0; end if (vull_V) t = [At: At: ant*At]; V = V(:,2:end); h = h(:,2:end); else 87