attempt to add lvad

.ipynb_checkpoints/app-checkpoint.py

@@ -495,7 +495,265 @@ def pvloop_simulator_plot_only(Rm, Ra, Emax, Emin, Vd, Tc, start_v):
     plt.title('Simulated PI-SSL LV Pressure Volume Loop', fontsize=16)
     return plot
 
-## Demo 
+#########################################
+# LVAD functions
+# RELU for diodes
+def r(u):
+    return max(u, 0.)
+
+def heart_ode0(y, t, Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd):
+    x1, x2, x3, x4, x5 = y #here y is a vector of 5 values (not functions), at time t, used for getting (dy/dt)(t)
+    P_lv = Plv(x1+Vd,Emax,Emin,t,Tc,Vd)
+    dydt = [r(x2-P_lv)/Rm-r(P_lv-x4)/Ra, (x3-x2)/(Rs*Cr)-r(x2-P_lv)/(Cr*Rm), (x2-x3)/(Rs*Cs)+x5/Cs, -x5/Ca+r(P_lv-x4)/(Ca*Ra), (x4-x3-Rc*x5)/Ls]
+    return dydt
+
+def getslope(y1, y2, y3, x1, x2, x3):
+  sum_x = x1 + x2 + x3
+  sum_y = y1 + y2 + y3
+  sum_xy = x1*y1 + x2*y2 + x3*y3
+  sum_xx = x1*x1 + x2*x2 + x3*x3
+  # calculate the coefficients of the least-squares line
+  n = 3
+  slope = (n*sum_xy - sum_x*sum_y) / (n*sum_xx - sum_x*sum_x)
+  return slope
+
+### ODE: for each t (here fixed), gives dy/dt as a function of y(t) at that t, so can be used for integrating the vector y over time
+#it is run for each t going from 0 to tmax
+def lvad_ode(y, t, Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew):
+    
+    #from simaan2008dynamical:
+    Ri = 0.0677
+    R0 = 0.0677
+    Rk = 0
+    x1bar = 1.
+    alpha = -3.5
+    Li = 0.0127
+    L0 = 0.0127
+    b0 = -0.296
+    b1 = -0.027
+    b2 = 9.9025e-7
+    
+    x1, x2, x3, x4, x5, x6, x7 = y #here y is a vector of 5 values (not functions), at time t, used for getting (dy/dt)(t)
+
+    P_lv = Plv(x1+Vd,Emax,Emin,t,Tc,Vd)
+    if (P_lv <= x1bar): Rk = alpha * (P_lv - x1bar)
+    Lstar = Li + L0 + b1
+    Lstar2 = -Li -L0 +b1
+    Rstar = Ri + + R0 + Rk + b0
+
+    dydt = [-x6 + r(x2-P_lv)/Rm-r(P_lv-x4)/Ra, (x3-x2)/(Rs*Cr)-r(x2-P_lv)/(Cr*Rm), (x2-x3)/(Rs*Cs)+x5/Cs, -x5/Ca+r(P_lv-x4)/(Ca*Ra) + x6/Ca, (x4-x3)/Ls-Rc*x5/Ls, -P_lv / Lstar2 + x4/Lstar2 + (Ri+R0+Rk-b0) / Lstar2 * x6 - b2 / Lstar2 * x7**2, ratew]
+
+    return dydt
+
+#returns pv loop and ef when there is no lvad:
+def f_nolvad(Tc, start_v, Emax, showpvloop):
+
+    N = 20
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    start_pla = float(start_v*Elastance(Emax, Emin, 0, Tc))
+    start_pao = 75.
+    start_pa = start_pao
+    start_qt = 0 #aortic flow is Q_T and is 0 at ED, also see Fig5 in simaan2008dynamical
+
+    y0 = [start_v, start_pla, start_pa, start_pao, start_qt]
+
+    t = np.linspace(0, Tc*N, int(60000*N)) #spaced numbers over interval (start, stop, number_of_steps), 60000 time instances for each heart cycle
+    #changed to 60000 for having integer positions for Tmax
+    #obtain 5D vector solution:
+    sol = odeint(heart_ode0, y0, t, args = (Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc,Vd)) #t: list of values
+
+    result_Vlv = np.array(sol[:, 0]) + Vd
+    result_Plv = np.array([Plv(v, Emax, Emin, xi, Tc, Vd) for xi,v in zip(t,sol[:, 0])])
+
+    #if showpvloop: plt.plot(result_Vlv[18*60000:20*60000], result_Plv[18*60000:20*60000], color='black', label='Without LVAD')
+
+    ved = sol[19*60000, 0] + Vd
+    ves = sol[200*int(60/Tc)+9000+19*60000, 0] + Vd
+    ef = (ved-ves)/ved * 100.
+    minv = min(result_Vlv[19*60000:20*60000-1])
+    minp = min(result_Plv[19*60000:20*60000-1])
+
+    result_pao = np.array(sol[:, 3])
+    pao_ed = min(result_pao[(N-1)*60000:N*60000-1])
+    pao_es = max(result_pao[(N-1)*60000:N*60000-1])
+
+    return ef, pao_ed, pao_es, ((ved - ves) * 60/Tc ) / 1000, sol[19*60000, 0], sol[19*60000, 1], sol[19*60000, 2], sol[19*60000, 3], sol[19*60000, 4], result_Vlv[18*60000:20*60000], result_Plv[18*60000:20*60000]
+
+#returns the w at which suction occurs: (i.e. for which the slope of the envelopes of x6 becomes negative)
+def get_suctionw(Tc, start_v, Emax, y00, y01, y02, y03, y04, w0, x60, ratew): #slope is slope0 for w
+
+    N = 70
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    start_pla = float(start_v*Elastance(Emax, Emin, 0, Tc))
+    start_pao = 75.
+    start_pa = start_pao
+    start_qt = 0 #aortic flow is Q_T and is 0 at ED, also see Fig5 in simaan2008dynamical
+
+    y0 = [start_v, start_pla, start_pa, start_pao, start_qt, x60, w0]
+    y0 = [y00, y01, y02, y03, y04, x60, w0]
+
+    ncycle = 20000
+    n = N * ncycle
+    sol = np.zeros((n, 7))
+    t = np.linspace(0., Tc * N, n)
+    for j in range(7):
+      sol[0][j] = y0[j]
+
+    result_Vlv = []
+    result_Plv = []
+    result_x6 = []
+    result_x7 = []
+    envx6 = []
+    timesenvx6 = []
+    slopes = []
+    ws = []
+
+    minx6 = 99999
+    tmin = 0
+    tlastupdate = 0
+    lastw = w0
+    update = 1
+
+    #solve the ODE step by step by adding dydt*dt:
+    for j in range(0, n-1):
+      #update y with dydt * dt
+      y = sol[j]
+      dydt = lvad_ode(y, t[j], Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew)
+      for k in range(7):
+        dydt[k] = dydt[k] * (t[j+1] - t[j])
+      sol[j+1] = sol[j] + dydt
+
+      #update the min of x6 in the current cylce. also keep the time at which the min is obtained (for getting the slope later)
+      if (minx6 > sol[j][5]):
+        minx6 = sol[j][5]
+        tmin = t[j]
+
+      #add minimum of x6 once each cycle ends: (works). then reset minx6 to 99999 for calculating again the minimum
+      if (j%ncycle==0 and j>1):
+        envx6.append(minx6)
+        timesenvx6.append(tmin)
+        minx6 = 99999
+
+      if (len(envx6)>=3):
+        slope = getslope(envx6[-1], envx6[-2], envx6[-3], timesenvx6[-1], timesenvx6[-2], timesenvx6[-3])
+        slopes.append(slope)
+        ws.append(y[6])
+
+    for i in range(n):
+      result_x6.append(sol[i, 5])
+      result_x7.append(sol[i, 6])
+
+    suction_w = 0
+    for i in range(2, len(slopes)):
+      if (slopes[i] < 0):
+        suction_w = ws[i-1]
+        break
+
+    return suction_w
+
+def f_lvad(Tc, start_v, Emax, c, slope, w0, x60, y00, y01, y02, y03, y04): #slope is slope0 for w
+
+    N = 70
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    y0 = [y00, y01, y02, y03, y04, x60, w0]
+
+    ncycle = 10000
+    n = N * ncycle
+    sol = np.zeros((n, 7))
+    t = np.linspace(0., Tc * N, n)
+    for j in range(7):
+      sol[0][j] = y0[j]
+
+    result_Vlv = []
+    result_Plv = []
+    result_x6 = []
+    result_x7 = []
+    envx6 = []
+    timesenvx6 = []
+
+    minx6 = 99999
+    tmin = 0
+    tlastupdate = 0
+    lastw = w0
+    update = 1
+    ratew = 0 #6000/60
+
+    #solve the ODE step by step by adding dydt*dt:
+    for j in range(0, n-1):
+      #update y with dydt * dt
+      y = sol[j]
+      dydt = lvad_ode(y, t[j], Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew)
+      for k in range(7):
+        dydt[k] = dydt[k] * (t[j+1] - t[j])
+      sol[j+1] = sol[j] + dydt
+
+      #update the min of x6 in the current cylce. also keep the time at which the min is obtained (for getting the slope later)
+      if (minx6 > sol[j][5]):
+        minx6 = sol[j][5]
+        tmin = t[j]
+
+      #add minimum of x6 once each cycle ends: (works). then reset minx6 to 99999 for calculating again the minimum
+      if (j%ncycle==0 and j>1):
+        envx6.append(minx6)
+        timesenvx6.append(tmin)
+        minx6 = 99999
+
+      #update w (if 0.005 s. have passed since the last update):
+      if (slope<0):
+        update = 0
+      if (t[j+1] - tlastupdate > 0.005 and slope>0 and update==1): #abs(slope)>0.0001
+        # if there are enough points of envelope: calculate slope:
+        if (len(envx6)>=3):
+          slope = getslope(envx6[-1], envx6[-2], envx6[-3], timesenvx6[-1], timesenvx6[-2], timesenvx6[-3])
+          sol[j+1][6] = lastw + c * slope
+        #otherwise: take arbitrary rate (see Fig. 16a in simaan2008dynamical)
+        else:
+          sol[j+1][6] = lastw + 0.005 * slope
+        #save w(k) (see formula (8) simaan2008dynamical) and the last time of update t[j] (will have to wait 0.005 s for next update of w)
+        tlastupdate = t[j+1]
+        lastw = sol[j+1][6]
+
+    #save functions and print MAP, CO:
+    map = 0
+    Pao = []
+
+    for i in range(n):
+      result_Vlv.append(sol[i, 0] + Vd)
+      result_Plv.append(Plv(sol[i, 0]+Vd, Emax, Emin, t[i], Tc, Vd))
+      result_x6.append(sol[i, 5])
+      result_x7.append(sol[i, 6])
+      Pao.append(sol[i, 3])
+
+    colors0=np.zeros((len(result_Vlv[65*ncycle:70*ncycle]), 3))
+    for col in colors0:
+      col[0]=41/255
+      col[1]=128/255
+      col[2]=205/255
+
+
+    #get co and ef:
+    ved = max(result_Vlv[50 * ncycle:52 * ncycle])
+    ves = min(result_Vlv[50 * ncycle:52 * ncycle])
+    #ves = result_Vlv[50 * ncycle + int(ncycle * 0.2 /Tc + 0.15 * ncycle)]
+    ef = (ved-ves)/ved*100
+    CO = ((ved - ves) * 60/Tc ) / 1000
+
+    #get MAP:
+    for i in range(n - 5*ncycle, n):
+      map += sol[i, 2]
+    map = map/(5*ncycle)
+
+    result_pao = np.array(sol[:, 3])
+    pao_ed = min(Pao[50 * ncycle:52 * ncycle])
+    pao_es = max(Pao[50 * ncycle:52 * ncycle])
+
+    return ef, pao_ed, pao_es, CO, map, result_Vlv[65*ncycle:70*ncycle], result_Plv[65*ncycle:70*ncycle]
+
+#############################
+## Demo functions
 
 def generate_example():
         # get random input
@@ -518,6 +776,53 @@ def generate_example():
         animated = "prediction.mp4"
         return video, animated, Rm, Ra, Emax, Emin, Vd, Tc, start_v
 
+
+def lvad_plots(Rm, Ra, Emax, Emin, Vd, Tc, start_v, gamma):
+
+    ncycle = 10000
+    
+    Rs = 1.
+    Rc = 0.0398
+    Ca= 0.08
+    Cs= 1.33
+    Cr= 4.4
+    Ls=0.0005
+    
+    #get values for periodic loops:
+    ef_nolvad, pao_ed, pao_es, co_nolvad, y00, y01, y02, y03, y04, Vlv0, Plv0 = f_nolvad(Tc, start_v, Emax, 0)
+    #pao_eds = [pao_ed]
+    #pao_ess = [pao_es]
+    
+    #get suction w: (make w go linearly from w0 to w0 + maxtime * 400, and find w at which suction occurs)
+    w0 = 5000.
+    ratew = 400.
+    x60 = 0.
+    suctionw = get_suctionw(Tc, start_v, Emax, y00, y01, y02, y03, y04, w0, x60, ratew)
+    
+    #gamma = 1.8
+    c = 0.065 #(in simaan2008dynamical: 0.67, but too fast -> 0.061 gives better shape)
+    slope0 = 100.
+    w0 = suctionw / gamma  #if doesn't work (x6 negative), change gamma down to 1.4 or up to 2.1
+    
+    #compute new pv loops and ef with lvad added:
+    new_ef, pao_ed, pao_es, CO, MAP, Vlvs, Plvs = f_lvad(Tc, start_v, Emax, c, slope0, w0, x60, y00, y01, y02, y03, y04)
+      
+    print("\nParameters: Tc, start_v, Emax:", Tc, start_v, Emax)
+    print("Suction speed:", suctionw)
+    print("EF before LVAD:", ef_nolvad)
+    print("CO before LVAD:", co_nolvad)
+    print("New EF after LVAD:", new_ef, "New CO:", CO, "New MAP:", MAP, "\n\n")
+    
+    plt.plot(Vlv0, Plv0, color='blue', label='No LVAD') #blue
+    plt.plot(Vlvs, Plvs, color=(78/255, 192/255, 44/255), label=f"LVAD, ω(0)= {round(w0,2)}r/min") #green
+    plt.xlabel('LV volume (ml)')
+    plt.ylabel('LV pressure (mmHg)')
+    plt.legend(loc='upper left', framealpha=1)
+    plt.ylim((0,220))
+    plt.xlim((0,250))
+    
+    return plt
+    
 title = "<h1 style='text-align: center; margin-bottom: 1rem'> Non-Invasive Medical Digital Twins using Physics-Informed Self-Supervised Learning </h1>"
 
 description = """
@@ -529,10 +834,16 @@ We demonstrate the ability of our model to predict left ventricular pressure-vol
 
 title2 = "<h3 style='text-align: center'> Physics-based model simulation</h3>"
 
-description2 = """\n \n
+description2 = """
+\n \n
 Our model uses a hydraulic analogy model of cardiac function from <a href='https://ieeexplore.ieee.org/document/4729737/keywords#keywords'>Simaan et al 2008</a>. Below you can input values of predicted parameters and output a simulated pressure-volume loop predicted from the <a href='https://ieeexplore.ieee.org/document/4729737/keywords#keywords'>Simaan et al 2008</a> model, which is an ordinary differential equation. Tune parameters and press 'Run simulation.'
 """
 
+description3 = """
+\n\n
+This model can incorporate a tunable left-ventricular assistance device (LVAD) for in-silico experimentation. Click to view the effect of adding an LVAD to the simulated PV loop.
+"""
+
 gr.Markdown(title)
 gr.Markdown(description)
 
@@ -573,13 +884,18 @@ with gr.Blocks() as demo:
                 sl5 = gr.Slider(4.0, 25.0, value= 4.0, label="Vd (ml)")
                 sl6 = gr.Slider(0.4, 1.7, value= 0.4, label="Tc (s)") 
                 sl7 = gr.Slider(0.0, 280.0, value= 140., label="start_v (ml)")
-    
+
+            gr.Markdown(description2)
+
+            LVAD_button = gr.Button("Add LVAD")
+
+            with gr.Row():
+                gamma = gr.Slider(1.0, 2.0, value= 1.4, label="Pump speed, ω(0)")
     
     generate_button.click(fn=generate_example, outputs = [video,plot,Rm,Ra,Emax,Emin,Vd,Tc,start_v])
     
-    
     simulation_button.click(fn=pvloop_simulator_plot_only, inputs = [sl1,sl2,sl3,sl4,sl5,sl6,sl7], outputs = [gr.Plot()])
     
-    
+    LVAD_button.click(fn=lvad_plots, inputs = [sl1,sl2,sl3,sl4,sl5,sl6,sl7,gamma], outputs = [gr.Plot()])
     
 demo.launch()

app.py

@@ -495,7 +495,265 @@ def pvloop_simulator_plot_only(Rm, Ra, Emax, Emin, Vd, Tc, start_v):
     plt.title('Simulated PI-SSL LV Pressure Volume Loop', fontsize=16)
     return plot
 
-## Demo 
+#########################################
+# LVAD functions
+# RELU for diodes
+def r(u):
+    return max(u, 0.)
+
+def heart_ode0(y, t, Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd):
+    x1, x2, x3, x4, x5 = y #here y is a vector of 5 values (not functions), at time t, used for getting (dy/dt)(t)
+    P_lv = Plv(x1+Vd,Emax,Emin,t,Tc,Vd)
+    dydt = [r(x2-P_lv)/Rm-r(P_lv-x4)/Ra, (x3-x2)/(Rs*Cr)-r(x2-P_lv)/(Cr*Rm), (x2-x3)/(Rs*Cs)+x5/Cs, -x5/Ca+r(P_lv-x4)/(Ca*Ra), (x4-x3-Rc*x5)/Ls]
+    return dydt
+
+def getslope(y1, y2, y3, x1, x2, x3):
+  sum_x = x1 + x2 + x3
+  sum_y = y1 + y2 + y3
+  sum_xy = x1*y1 + x2*y2 + x3*y3
+  sum_xx = x1*x1 + x2*x2 + x3*x3
+  # calculate the coefficients of the least-squares line
+  n = 3
+  slope = (n*sum_xy - sum_x*sum_y) / (n*sum_xx - sum_x*sum_x)
+  return slope
+
+### ODE: for each t (here fixed), gives dy/dt as a function of y(t) at that t, so can be used for integrating the vector y over time
+#it is run for each t going from 0 to tmax
+def lvad_ode(y, t, Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew):
+    
+    #from simaan2008dynamical:
+    Ri = 0.0677
+    R0 = 0.0677
+    Rk = 0
+    x1bar = 1.
+    alpha = -3.5
+    Li = 0.0127
+    L0 = 0.0127
+    b0 = -0.296
+    b1 = -0.027
+    b2 = 9.9025e-7
+    
+    x1, x2, x3, x4, x5, x6, x7 = y #here y is a vector of 5 values (not functions), at time t, used for getting (dy/dt)(t)
+
+    P_lv = Plv(x1+Vd,Emax,Emin,t,Tc,Vd)
+    if (P_lv <= x1bar): Rk = alpha * (P_lv - x1bar)
+    Lstar = Li + L0 + b1
+    Lstar2 = -Li -L0 +b1
+    Rstar = Ri + + R0 + Rk + b0
+
+    dydt = [-x6 + r(x2-P_lv)/Rm-r(P_lv-x4)/Ra, (x3-x2)/(Rs*Cr)-r(x2-P_lv)/(Cr*Rm), (x2-x3)/(Rs*Cs)+x5/Cs, -x5/Ca+r(P_lv-x4)/(Ca*Ra) + x6/Ca, (x4-x3)/Ls-Rc*x5/Ls, -P_lv / Lstar2 + x4/Lstar2 + (Ri+R0+Rk-b0) / Lstar2 * x6 - b2 / Lstar2 * x7**2, ratew]
+
+    return dydt
+
+#returns pv loop and ef when there is no lvad:
+def f_nolvad(Tc, start_v, Emax, showpvloop):
+
+    N = 20
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    start_pla = float(start_v*Elastance(Emax, Emin, 0, Tc))
+    start_pao = 75.
+    start_pa = start_pao
+    start_qt = 0 #aortic flow is Q_T and is 0 at ED, also see Fig5 in simaan2008dynamical
+
+    y0 = [start_v, start_pla, start_pa, start_pao, start_qt]
+
+    t = np.linspace(0, Tc*N, int(60000*N)) #spaced numbers over interval (start, stop, number_of_steps), 60000 time instances for each heart cycle
+    #changed to 60000 for having integer positions for Tmax
+    #obtain 5D vector solution:
+    sol = odeint(heart_ode0, y0, t, args = (Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc,Vd)) #t: list of values
+
+    result_Vlv = np.array(sol[:, 0]) + Vd
+    result_Plv = np.array([Plv(v, Emax, Emin, xi, Tc, Vd) for xi,v in zip(t,sol[:, 0])])
+
+    #if showpvloop: plt.plot(result_Vlv[18*60000:20*60000], result_Plv[18*60000:20*60000], color='black', label='Without LVAD')
+
+    ved = sol[19*60000, 0] + Vd
+    ves = sol[200*int(60/Tc)+9000+19*60000, 0] + Vd
+    ef = (ved-ves)/ved * 100.
+    minv = min(result_Vlv[19*60000:20*60000-1])
+    minp = min(result_Plv[19*60000:20*60000-1])
+
+    result_pao = np.array(sol[:, 3])
+    pao_ed = min(result_pao[(N-1)*60000:N*60000-1])
+    pao_es = max(result_pao[(N-1)*60000:N*60000-1])
+
+    return ef, pao_ed, pao_es, ((ved - ves) * 60/Tc ) / 1000, sol[19*60000, 0], sol[19*60000, 1], sol[19*60000, 2], sol[19*60000, 3], sol[19*60000, 4], result_Vlv[18*60000:20*60000], result_Plv[18*60000:20*60000]
+
+#returns the w at which suction occurs: (i.e. for which the slope of the envelopes of x6 becomes negative)
+def get_suctionw(Tc, start_v, Emax, y00, y01, y02, y03, y04, w0, x60, ratew): #slope is slope0 for w
+
+    N = 70
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    start_pla = float(start_v*Elastance(Emax, Emin, 0, Tc))
+    start_pao = 75.
+    start_pa = start_pao
+    start_qt = 0 #aortic flow is Q_T and is 0 at ED, also see Fig5 in simaan2008dynamical
+
+    y0 = [start_v, start_pla, start_pa, start_pao, start_qt, x60, w0]
+    y0 = [y00, y01, y02, y03, y04, x60, w0]
+
+    ncycle = 20000
+    n = N * ncycle
+    sol = np.zeros((n, 7))
+    t = np.linspace(0., Tc * N, n)
+    for j in range(7):
+      sol[0][j] = y0[j]
+
+    result_Vlv = []
+    result_Plv = []
+    result_x6 = []
+    result_x7 = []
+    envx6 = []
+    timesenvx6 = []
+    slopes = []
+    ws = []
+
+    minx6 = 99999
+    tmin = 0
+    tlastupdate = 0
+    lastw = w0
+    update = 1
+
+    #solve the ODE step by step by adding dydt*dt:
+    for j in range(0, n-1):
+      #update y with dydt * dt
+      y = sol[j]
+      dydt = lvad_ode(y, t[j], Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew)
+      for k in range(7):
+        dydt[k] = dydt[k] * (t[j+1] - t[j])
+      sol[j+1] = sol[j] + dydt
+
+      #update the min of x6 in the current cylce. also keep the time at which the min is obtained (for getting the slope later)
+      if (minx6 > sol[j][5]):
+        minx6 = sol[j][5]
+        tmin = t[j]
+
+      #add minimum of x6 once each cycle ends: (works). then reset minx6 to 99999 for calculating again the minimum
+      if (j%ncycle==0 and j>1):
+        envx6.append(minx6)
+        timesenvx6.append(tmin)
+        minx6 = 99999
+
+      if (len(envx6)>=3):
+        slope = getslope(envx6[-1], envx6[-2], envx6[-3], timesenvx6[-1], timesenvx6[-2], timesenvx6[-3])
+        slopes.append(slope)
+        ws.append(y[6])
+
+    for i in range(n):
+      result_x6.append(sol[i, 5])
+      result_x7.append(sol[i, 6])
+
+    suction_w = 0
+    for i in range(2, len(slopes)):
+      if (slopes[i] < 0):
+        suction_w = ws[i-1]
+        break
+
+    return suction_w
+
+def f_lvad(Tc, start_v, Emax, c, slope, w0, x60, y00, y01, y02, y03, y04): #slope is slope0 for w
+
+    N = 70
+    global Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emin, Vd
+
+    y0 = [y00, y01, y02, y03, y04, x60, w0]
+
+    ncycle = 10000
+    n = N * ncycle
+    sol = np.zeros((n, 7))
+    t = np.linspace(0., Tc * N, n)
+    for j in range(7):
+      sol[0][j] = y0[j]
+
+    result_Vlv = []
+    result_Plv = []
+    result_x6 = []
+    result_x7 = []
+    envx6 = []
+    timesenvx6 = []
+
+    minx6 = 99999
+    tmin = 0
+    tlastupdate = 0
+    lastw = w0
+    update = 1
+    ratew = 0 #6000/60
+
+    #solve the ODE step by step by adding dydt*dt:
+    for j in range(0, n-1):
+      #update y with dydt * dt
+      y = sol[j]
+      dydt = lvad_ode(y, t[j], Rs, Rm, Ra, Rc, Ca, Cs, Cr, Ls, Emax, Emin, Tc, Vd, ratew)
+      for k in range(7):
+        dydt[k] = dydt[k] * (t[j+1] - t[j])
+      sol[j+1] = sol[j] + dydt
+
+      #update the min of x6 in the current cylce. also keep the time at which the min is obtained (for getting the slope later)
+      if (minx6 > sol[j][5]):
+        minx6 = sol[j][5]
+        tmin = t[j]
+
+      #add minimum of x6 once each cycle ends: (works). then reset minx6 to 99999 for calculating again the minimum
+      if (j%ncycle==0 and j>1):
+        envx6.append(minx6)
+        timesenvx6.append(tmin)
+        minx6 = 99999
+
+      #update w (if 0.005 s. have passed since the last update):
+      if (slope<0):
+        update = 0
+      if (t[j+1] - tlastupdate > 0.005 and slope>0 and update==1): #abs(slope)>0.0001
+        # if there are enough points of envelope: calculate slope:
+        if (len(envx6)>=3):
+          slope = getslope(envx6[-1], envx6[-2], envx6[-3], timesenvx6[-1], timesenvx6[-2], timesenvx6[-3])
+          sol[j+1][6] = lastw + c * slope
+        #otherwise: take arbitrary rate (see Fig. 16a in simaan2008dynamical)
+        else:
+          sol[j+1][6] = lastw + 0.005 * slope
+        #save w(k) (see formula (8) simaan2008dynamical) and the last time of update t[j] (will have to wait 0.005 s for next update of w)
+        tlastupdate = t[j+1]
+        lastw = sol[j+1][6]
+
+    #save functions and print MAP, CO:
+    map = 0
+    Pao = []
+
+    for i in range(n):
+      result_Vlv.append(sol[i, 0] + Vd)
+      result_Plv.append(Plv(sol[i, 0]+Vd, Emax, Emin, t[i], Tc, Vd))
+      result_x6.append(sol[i, 5])
+      result_x7.append(sol[i, 6])
+      Pao.append(sol[i, 3])
+
+    colors0=np.zeros((len(result_Vlv[65*ncycle:70*ncycle]), 3))
+    for col in colors0:
+      col[0]=41/255
+      col[1]=128/255
+      col[2]=205/255
+
+
+    #get co and ef:
+    ved = max(result_Vlv[50 * ncycle:52 * ncycle])
+    ves = min(result_Vlv[50 * ncycle:52 * ncycle])
+    #ves = result_Vlv[50 * ncycle + int(ncycle * 0.2 /Tc + 0.15 * ncycle)]
+    ef = (ved-ves)/ved*100
+    CO = ((ved - ves) * 60/Tc ) / 1000
+
+    #get MAP:
+    for i in range(n - 5*ncycle, n):
+      map += sol[i, 2]
+    map = map/(5*ncycle)
+
+    result_pao = np.array(sol[:, 3])
+    pao_ed = min(Pao[50 * ncycle:52 * ncycle])
+    pao_es = max(Pao[50 * ncycle:52 * ncycle])
+
+    return ef, pao_ed, pao_es, CO, map, result_Vlv[65*ncycle:70*ncycle], result_Plv[65*ncycle:70*ncycle]
+
+#############################
+## Demo functions
 
 def generate_example():
         # get random input
@@ -518,6 +776,53 @@ def generate_example():
         animated = "prediction.mp4"
         return video, animated, Rm, Ra, Emax, Emin, Vd, Tc, start_v
 
+
+def lvad_plots(Rm, Ra, Emax, Emin, Vd, Tc, start_v, gamma):
+
+    ncycle = 10000
+    
+    Rs = 1.
+    Rc = 0.0398
+    Ca= 0.08
+    Cs= 1.33
+    Cr= 4.4
+    Ls=0.0005
+    
+    #get values for periodic loops:
+    ef_nolvad, pao_ed, pao_es, co_nolvad, y00, y01, y02, y03, y04, Vlv0, Plv0 = f_nolvad(Tc, start_v, Emax, 0)
+    #pao_eds = [pao_ed]
+    #pao_ess = [pao_es]
+    
+    #get suction w: (make w go linearly from w0 to w0 + maxtime * 400, and find w at which suction occurs)
+    w0 = 5000.
+    ratew = 400.
+    x60 = 0.
+    suctionw = get_suctionw(Tc, start_v, Emax, y00, y01, y02, y03, y04, w0, x60, ratew)
+    
+    #gamma = 1.8
+    c = 0.065 #(in simaan2008dynamical: 0.67, but too fast -> 0.061 gives better shape)
+    slope0 = 100.
+    w0 = suctionw / gamma  #if doesn't work (x6 negative), change gamma down to 1.4 or up to 2.1
+    
+    #compute new pv loops and ef with lvad added:
+    new_ef, pao_ed, pao_es, CO, MAP, Vlvs, Plvs = f_lvad(Tc, start_v, Emax, c, slope0, w0, x60, y00, y01, y02, y03, y04)
+      
+    print("\nParameters: Tc, start_v, Emax:", Tc, start_v, Emax)
+    print("Suction speed:", suctionw)
+    print("EF before LVAD:", ef_nolvad)
+    print("CO before LVAD:", co_nolvad)
+    print("New EF after LVAD:", new_ef, "New CO:", CO, "New MAP:", MAP, "\n\n")
+    
+    plt.plot(Vlv0, Plv0, color='blue', label='No LVAD') #blue
+    plt.plot(Vlvs, Plvs, color=(78/255, 192/255, 44/255), label=f"LVAD, ω(0)= {round(w0,2)}r/min") #green
+    plt.xlabel('LV volume (ml)')
+    plt.ylabel('LV pressure (mmHg)')
+    plt.legend(loc='upper left', framealpha=1)
+    plt.ylim((0,220))
+    plt.xlim((0,250))
+    
+    return plt
+    
 title = "<h1 style='text-align: center; margin-bottom: 1rem'> Non-Invasive Medical Digital Twins using Physics-Informed Self-Supervised Learning </h1>"
 
 description = """
@@ -529,10 +834,16 @@ We demonstrate the ability of our model to predict left ventricular pressure-vol
 
 title2 = "<h3 style='text-align: center'> Physics-based model simulation</h3>"
 
-description2 = """\n \n
+description2 = """
+\n \n
 Our model uses a hydraulic analogy model of cardiac function from <a href='https://ieeexplore.ieee.org/document/4729737/keywords#keywords'>Simaan et al 2008</a>. Below you can input values of predicted parameters and output a simulated pressure-volume loop predicted from the <a href='https://ieeexplore.ieee.org/document/4729737/keywords#keywords'>Simaan et al 2008</a> model, which is an ordinary differential equation. Tune parameters and press 'Run simulation.'
 """
 
+description3 = """
+\n\n
+This model can incorporate a tunable left-ventricular assistance device (LVAD) for in-silico experimentation. Click to view the effect of adding an LVAD to the simulated PV loop.
+"""
+
 gr.Markdown(title)
 gr.Markdown(description)
 
@@ -573,12 +884,18 @@ with gr.Blocks() as demo:
                 sl5 = gr.Slider(4.0, 25.0, value= 4.0, label="Vd (ml)")
                 sl6 = gr.Slider(0.4, 1.7, value= 0.4, label="Tc (s)") 
                 sl7 = gr.Slider(0.0, 280.0, value= 140., label="start_v (ml)")
-    
+
+            gr.Markdown(description2)
+
+            LVAD_button = gr.Button("Add LVAD")
+
+            with gr.Row():
+                gamma = gr.Slider(1.0, 2.0, value= 1.4, label="Pump speed, ω(0)")
     
     generate_button.click(fn=generate_example, outputs = [video,plot,Rm,Ra,Emax,Emin,Vd,Tc,start_v])
     
     simulation_button.click(fn=pvloop_simulator_plot_only, inputs = [sl1,sl2,sl3,sl4,sl5,sl6,sl7], outputs = [gr.Plot()])
     
-    
+    LVAD_button.click(fn=lvad_plots, inputs = [sl1,sl2,sl3,sl4,sl5,sl6,sl7,gamma], outputs = [gr.Plot()])
     
 demo.launch()