Home | Scholars at Harvard



Peter Chang2016 Summer Research SummarySummaryIn the deep ocean, the tunicate, more commonly known as a sea squirt, exhibits peculiar behavior with its heart. While all other animals have hearts that pump blood in one direction, the sea squirt is able to pump blood in both directions. Neither the biological incentive nor the physiological mechanism has been well studied, however there are potential medical benefits to understanding this phenomena.The goal of this project was to determine potential mathematical methods that govern the reversals within the heart, which would be tested through simulations. The two tested methods that both produced the expected results were a random variation method and a controlled shifting method. These results, as well as the mechanisms behind them, will be presented in this summary.Basic Physiological EquationsIn order for a heart to pump blood, a pacemaker is required at the end of the heart fibers. This pacemaker creates electric jolts at a certain interval in order to send waves throughout the entire fiber. The heart of a sea squirt may be modeled as having two pacemakers, one at either end of the heart fiber (Laura Miller), which allows for blood to flow in both directions. There are two main components that govern the propagation of waves within the heart, the first is the processes of individual heart cells and the second is diffusion between adjacent heart cells.When dealing with individual heart cells, there are two differential equations that govern how electric potential is stored (Mitchell and Schaeffer). The first is the primary equation for voltage:dvdt=inward current-outward current=htinv21-v-vtoutWhere tin and tout are physiological constants, v is voltage inside the cell, and h is a constant between 0 and 1 that represents how open or closed the cell door which introduces voltage from outside the fiber. The variable h is governed by the ordinary differential equation:dhdt=-htclose if v>vcrit=1-htopen if v<vcritWhere topen and tclose are once again physiological constants and vcrit is a certain voltage level that determines how the cell will behave.The second important component is the diffusion between cells. The diffusion equation for voltage with respect to time is:dudt=kd2udx2Where k is another physiological constant. Diffusion will be a factor of the cell itself, and the cells immediately adjacent to it. We can estimate the change by diffusion with the following equation:Kuk-1m-2ukm+uk-1m(?x)2Where K is the diffusion coefficient (.001 for most cases), and ukm is the estimated voltage at time m and cell k. There are boundary conditions on the cells at the beginning and ends of the fiber. By assuming the conditions that no voltage passes through the ends, i.e. dudx0,t=0 and dudxL,t=0, these will have an increased effect from the single adjacent cell:K2u1m-2u0m(?x)2K-2uNm+2uN-1m(?x)2With N representing the number of cells within the heart fiber.Approximation MethodsSince we are dealing with differential equations, we had to use an approximation method in order to compute the necessary values. Our initial goal was to implement a Crank-Nicolson solver for this, however there were some unexpected outputs with the order of convergence of the error estimations. Overall, the Crank-Nicolson solver seemed to be accurate, however we decided not to use it as a safety measure.Instead, we implemented a forward Euler solver. This method requires a smaller timestep during the calculations, which meant that it was much slower than the Crank-Nicolson solver. However, this increased computation time was necessary in order to ensure the stability and robustness of the code.In its most general form, the forward Euler approximation is given an initial condition and calculates each consecutive value with:yn+1=yn+?t*f(tn, yn)If we apply this to the equations that we have already outlined above, and combine the individual cell and diffusion equations together, we get the overall equations for computing the voltage, v, and door constant, h, at a location k and time m+1:vkm+1=vkm+?thv21-vtin-vtout+Kvk+1m-2vkm+vk-1m (?x)2*hkm+1=hkm+?t1-htopen if v≤ vcrit -htcloseif v>vcrit*This part will change slightly for cells at the end of the fiber to account for boundary conditions.Thus we have two equations that govern the propagation of the waves through the fiber. In order to get this as a full model, however, there also needs to be initial conditions. In a completely rested heart fiber, the door constant, h, at every cell would be set to 1 and the voltages would all be 0. We can replicate the pacemakers by “jolting” the cells at the end of the fiber to a certain voltage. Even just by instantly raising the voltage of the first four cells (out of several hundred), we created a wave which would propagate through the rest of the cells, with a sharp wall leading the way and the rest of the cells “resting” down to no voltage.Figure 1: single wave propagatingCoding ComponentThese equations were coded into Python in order to generate simulations of the phenomena. For this part, there were two arrays, one representing the voltage of each individual cell, and one representing the door constant h for each individual cell. Using the equations above, the arrays could continuously be updated and the data stored in it could be exported as a graph. These graphs were then stitched together using ffmpeg to generate full movies.Saving memory, time, and CPU usage were all major concerns throughout this process. An initial code which did not erase old arrays as they were updated was extremely slow and resulted in memory overloads. By exporting the graphs only every couple hundred updates and deleting the old voltage, the code was able to use a fraction of that memory and CPU. Python proved to be very clunky in this project, as equivalent code in C++ tended to run orders of magnitude faster. Hypothesis TestingSeveral physiological constants were needed to complete the equations. Constants like K = .001 and vcrit = .13 are constants that are widely accepted in the modeling of heart cells. The other constants, tin = .1, tout = 2.4, topen = 150, and tclose = 180, were all set to match human heart fiber constants. These were slightly adjusted and the different results were recorded throughout the testing process to account for tunicate heart cells being slightly different from human heart cells.In addition, it was important to set the values for length, number of cells, and computational constants. We used a length of 3 centimeters (a common length for tunicate heart fibers) with 300 individual cells, making each cell .1 millimeters wide. All time was measured in milliseconds and computations were done every .01ms. These computational constants were chosen so the CFL condition of K?t?x2<12 was easily satisfied, which meant that our computations would be stable under the forward Euler approximation.When two waves coming from opposite ends of the fiber meet, they will stop each other at the location they meet, and slowly die down together. This is a feature that differentiates the waves in the heart versus waves in most other medium and is a result of the conditions that are present between each cell. This means that if waves were propagating at a faster rate from one side of the fiber, that side would completely dominate all blood flow within the heart.To first generate reversals within the fibers, each pacemaker at the ends were set to different beating rates. The left side was set to pulse every 500ms and after each pulse, the interval would be lowered by 10ms (so the second pulse would come after another 490ms). The right side was set to pulse every 350ms and similarly, the interval would increase after 10ms. Once these two pacemakers have completely switched their paces (i.e. left side at 350ms and right side at 500ms), they would reverse direction.As would clearly be expected, this did generate reversals, however they were not as clear as one might think. The side with the shorter interval did indeed dominate the direction of the blood flow, however there were occasional waves from the opposite side that would sneak through. In general, however, reversals would only occur once the intervals of the two pacemakers had changed dominance.The second method tested whether two pacemakers at the same pace could generate reversals through only natural variation in the beating. For this method, the two pacemakers were set to beat on average at the same pace, every 400ms, however the individual intervals fluctuated with a standard deviation of 20ms. This method also generated reversals and these reversals were much cleaner than that of the first method. That is, there were no random waves which would sneak through and temporarily reverse the direction of the blood flow. In this random fluctuation method, the reversals only occurred after long periods of time had elapsed and occurred very rapidly, which is consistent with tunicate behavior in nature. It is important to note that this method generated fast reversals when the fiber was of length 3cm like in our simulations. When testing on a 10cm fiber, the reversals took much longer since the random fluctuations would take much longer to generate a new dominant side. This, however, is not much of an issue since the 3cm length is much more accurate for a tunicate.ConclusionThrough the testing done, it seems clear that both of these methods are viable options for the reversal of blood flow in tunicate hearts. When viewing the simulations, it is easy to see the reversal of blood flow happening. While both of them did produce the effect, the method that used pacemakers with equal interval but with some random variation (which is typical in any real world heart) produced a cleaner and more realistic result.For further testing, it would be nice to convert the code into C++ or some other language that runs faster computations that Python. It would also have been good to try to fix the Crank-Nicolson solver, which also would have decreased the computation time. These fixes would have allowed for more simulations to be generated so that more hypotheses could have been tested. ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download