-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotter.py
287 lines (239 loc) · 10.5 KB
/
plotter.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy as np
from fenics import *
# set font & text parameters
font = {"family": "serif", "weight": "bold", "size": 13}
plt.rc("font", **font)
plt.rc("text", usetex=True)
mpl.rcParams["image.cmap"] = "jet"
class Plotter:
def __init__(self, res_2D, T, dt, f, f_emi, f_hyper, f_emi_hyper):
# plotting parameters
self.res_2D = res_2D # mesh resolution
self.T = T # end time
self.dt = 1.0e-3
# create time series
self.time = 1.0e3 * np.arange(dt, T, dt)
# files containing solutions
self.f = f # KNP-EMI normal activity
self.f_hyper = f_hyper # KNP-EMI hyperactivity
self.f_emi = f_emi # EMI normal activity
self.f_emi_hyper = f_emi_hyper # EMI normal activity
return
def get_plottable_ECS_function(self, h5_fname, n, i, scale=1.0):
""" Get plottable function of extracellular concentration or potential """
mesh = Mesh()
subdomains = MeshFunction("size_t", mesh, 2)
surfaces = MeshFunction("size_t", mesh, 1)
hdf5 = HDF5File(MPI.comm_world, h5_fname, "r")
hdf5.read(mesh, "/mesh", False)
mesh.coordinates()[:] *= 1e6
hdf5.read(subdomains, "/subdomains")
hdf5.read(surfaces, "/surfaces")
exterior_mesh = SubMesh(mesh, subdomains, 0)
P1 = FiniteElement("CG", triangle, 1)
R = FiniteElement("R", triangle, 0)
# EMI
if i is None:
We = FunctionSpace(exterior_mesh, MixedElement([P1] + [R]))
# KNP-EMI
else:
We = FunctionSpace(exterior_mesh, MixedElement(4 * [P1] + [R]))
Ve = FunctionSpace(exterior_mesh, P1)
ue = Function(We)
fe = Function(Ve)
# EMI
if i is None:
hdf5.read(ue, "/exterior_solution/vector_" + str(n))
assign(fe, ue.sub(0))
# KNP-EMI
else:
hdf5.read(ue, "/exterior_solution/vector_" + str(n))
assign(fe, ue.sub(i))
# scale (e.g. from V to mV)
fe.vector()[:] = scale * fe.vector().get_local()
return fe
def get_time_series_ECS(self, dt, T, fname, x, y, z, EMI=False):
""" Return list of values in given point (x, y, z) over time """
# read data file
hdf5file = HDF5File(MPI.comm_world, fname, "r")
# create mesh
mesh = Mesh()
subdomains = MeshFunction("size_t", mesh, 2)
surfaces = MeshFunction("size_t", mesh, 1)
hdf5file.read(mesh, "/mesh", False)
mesh.coordinates()[:] *= 1e6
hdf5file.read(subdomains, "/subdomains")
hdf5file.read(surfaces, "/surfaces")
exterior_mesh = SubMesh(mesh, subdomains, 0)
# define function spaces
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement('R', mesh.ufl_cell(), 0)
# define functions
Ve = FunctionSpace(exterior_mesh, P1)
f_phi_e = Function(Ve)
f_Na_e = Function(Ve)
f_K_e = Function(Ve)
# list of values at point over time
Na_e = np.empty(int(T / dt) - 1)
K_e = np.empty(int(T / dt) - 1)
phi_e = np.empty(int(T / dt) - 1)
for n in range(1, int(T / dt)):
if EMI:
# create function space
We = FunctionSpace(exterior_mesh, MixedElement([P1] + [R]))
ue = Function(We)
# read file
hdf5file.read(ue, "/exterior_solution/vector_" + str(n))
# assign potential
assign(f_phi_e, ue.sub(0))
# append potential
phi_e[n - 1] = 1.0e3 * f_phi_e(x, y, z) # 3D
else:
# create function space
We = FunctionSpace(exterior_mesh, MixedElement(4 * [P1] + [R]))
ue = Function(We)
# read file
hdf5file.read(ue, "/exterior_solution/vector_" + str(n))
# assign potentials and concentrations
assign(f_Na_e, ue.sub(0)) # ECS Na concentrations
assign(f_K_e, ue.sub(1)) # ECS K concentrations
assign(f_phi_e, ue.sub(3)) # ECS potential
# append potentials and concentrations
phi_e[n - 1] = 1.0e3 * f_phi_e(x, y, z)
Na_e[n - 1] = f_Na_e(x, y, z)
K_e[n - 1] = f_K_e(x, y, z)
return Na_e, K_e, phi_e
def get_time_series_gamma(self, dt, T, fname, x, y, z):
""" Return list of values in given point (x, y, z) over time """
# read data file
hdf5file = HDF5File(MPI.comm_world, fname, "r")
# membrane potential
gamma_mesh = Mesh()
hdf5file.read(gamma_mesh, "/gamma_mesh", False)
P1 = FiniteElement("P", gamma_mesh.ufl_cell(), 1)
Vg = FunctionSpace(gamma_mesh, P1)
gamma_mesh.coordinates()[:] *= 1e6
f_phi_M = Function(Vg)
phi_M = []
for n in range(1, int(T / dt)):
# read file
hdf5file.read(f_phi_M, "/membrane_potential/vector_" + str(n))
# membrane potential
phi_M.append(1.0e3 * f_phi_M(x, y, z))
return phi_M
def make_figures(self):
""" Create plots of potentials and ion concentrations """
# set time parameters
dt = self.dt
T = self.T
time = self.time
# point at membrane
x_M = 100; y_M = 0.5; z_M = 0.5
# point in ECS (1 um above axon)
x_E = 100; y_E = 0.5 + 0.05; z_E = 0.5
# get time series - membrane potential
phi_M = self.get_time_series_gamma(dt, T, self.f, x_M, y_M, z_M)
phi_M_emi = self.get_time_series_gamma(dt, T, self.f_emi, x_M, y_M, z_M)
# get time series - ECS ion concentrations
Na, K, phi_E = self.get_time_series_ECS(dt, T, self.f, x_E, y_E, z_E)
_, _, phi_E_emi = self.get_time_series_ECS(dt, T, self.f_emi, x_E, y_E, z_E, EMI=True)
# get time series - membrane potential
phi_M_hyper = self.get_time_series_gamma(dt, T, self.f_hyper, x_M, y_M, z_M)
phi_M_emi_hyper = self.get_time_series_gamma(dt, T, self.f_emi_hyper, x_M, y_M, z_M)
# get time series - ion concentrations
Na_hyper, K_hyper, phi_E_hyper = self.get_time_series_ECS(dt, T, self.f_hyper, x_E, y_E, z_E)
_, _, phi_E_emi_hyper = self.get_time_series_ECS(dt, T, self.f_emi_hyper, x_E, y_E, z_E, EMI=True)
# create figure for potential
fig = plt.figure(figsize=(12 / 1.5, 25 / 3.0))
gs = GridSpec(nrows=2, ncols=2, width_ratios=[1, 1], height_ratios=[1, 1])
gs.update(left=0.1, right=0.95, bottom=0.1, top=0.93, wspace=0.3, hspace=0.3)
# membrane potential during normal activity
ax1 = fig.add_subplot(gs[0, 0])
plt.title(r"Membrane potential")
plt.ylabel(r"$v$ (mV)")
plt.xlabel(r"time (ms)")
plt.yticks([-80, -60, -40, -20, 0, 20])
plt.ylim(-100, 40)
plt.plot(time, phi_M, linewidth=3, label="KNP-EMI")
plt.plot(time, phi_M_emi, "--", linewidth=3, label="EMI")
plt.legend()
# membrane potential during hyperactivity
ax2 = fig.add_subplot(gs[1, 0])
plt.title(r"Membrane potential")
plt.ylabel(r"$v$ (mV)")
plt.xlabel(r"time (ms)")
plt.yticks([-80, -60, -40, -20, 0, 20])
plt.ylim(-100, 40)
plt.plot(time, phi_M_hyper, linewidth=3, label="KNP-EMI")
plt.plot(time, phi_M_emi_hyper, "--", linewidth=3, label="EMI")
# ECS potential during normal activity
ax3 = fig.add_subplot(gs[0, 1])
plt.title(r"ECS potential")
plt.ylabel(r"$u_e$ (mV)")
plt.xlabel(r"time (ms)")
plt.yticks([-0.8, -0.4, 0, 0.4, 0.8])
plt.ylim(-1.0, 1.0)
plt.plot(time, phi_E, linewidth=3, label="KNP-EMI")
plt.plot(time, phi_E_emi, "--", linewidth=3, label="EMI")
# ECS potential during hyperactivity
ax4 = fig.add_subplot(gs[1, 1])
plt.title(r"ECS potential")
plt.ylabel(r"$u_e$ (mV)")
plt.xlabel(r"time (ms)")
plt.yticks([-0.8, -0.4, 0, 0.4, 0.8])
plt.ylim(-1.0, 1.0)
plt.plot(time, phi_E_hyper, linewidth=3, label="KNP-EMI")
plt.plot(time, phi_E_emi_hyper, "--", linewidth=3, label="EMI")
# add numbering for the subplots (A, B, C etc)
letters = [r"\textbf{A}", r"\textbf{C}", r"\textbf{B}", r"\textbf{D}"]
for n, ax in enumerate([ax1, ax2, ax3, ax4]):
ax.text(-0.1, 1.1, letters[n], transform=ax.transAxes)
# save figure
plt.savefig("potentials.svg")
# create figure for ECS concentrations
fig = plt.figure(figsize=(12 / 1.5, 25 / 3.0))
gs = GridSpec(nrows=2, ncols=2, width_ratios=[1, 1], height_ratios=[1, 1])
gs.update(left=0.1, right=0.95, bottom=0.1, top=0.93, wspace=0.3, hspace=0.3)
# ECS Na concentration during normal activity
ax1 = fig.add_subplot(gs[0, 0])
plt.title(r"ECS Na$^+$ concentration")
plt.ylabel(r"$c^{Na}_e$(mM)")
plt.xlabel(r"time (ms)")
plt.yticks([119.4, 119.6, 119.8, 120])
plt.ylim([119.3, 120.1])
plt.plot(time, Na, linewidth=3)
# ECS Na concentration during hyperactivity
ax2 = fig.add_subplot(gs[1, 0])
plt.title(r"ECS Na$^+$ concentration")
plt.ylabel(r"$c^{Na}_e$(mM)")
plt.xlabel(r"time (ms)")
plt.yticks([119.4, 119.6, 119.8, 120])
plt.ylim([119.3, 120.1])
plt.plot(time, Na_hyper, linewidth=3)
# ECS K concentration during normal activity
ax3 = fig.add_subplot(gs[0, 1])
plt.title(r"ECS K$^+$ concentration")
plt.ylabel(r"$c^{K}_e$(mM)")
plt.xlabel(r"time (ms)")
plt.yticks([4.0, 4.2, 4.4, 4.6, 4.8])
plt.ylim([3.96, 4.84])
plt.plot(time, K, linewidth=3)
# ECS K concentration during hyperactivity
ax4 = fig.add_subplot(gs[1, 1])
plt.title(r"ECS K$^+$ concentration")
plt.ylabel(r"$c^{K}_e$(mM)")
plt.xlabel(r"time (ms)")
plt.yticks([4.0, 4.2, 4.4, 4.6, 4.8])
plt.ylim([3.96, 4.84])
plt.plot(time, K_hyper, linewidth=3)
# add numbering for the subplots (A, B, C etc)
letters = [r"\textbf{A}", r"\textbf{C}", r"\textbf{B}", r"\textbf{D}"]
for n, ax in enumerate([ax1, ax2, ax3, ax4]):
ax.text(-0.1, 1.1, letters[n], transform=ax.transAxes)
# save figure
plt.savefig("concentrations.svg")
return