Coverage for physioblocks / library / model_components / active_law.py: 100%

1# SPDX-FileCopyrightText: Copyright INRIA 

2# 

3# SPDX-License-Identifier: LGPL-3.0-only 

4# 

5# Copyright INRIA 

6# 

7# This file is part of PhysioBlocks, a library mostly developed by the 

8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA. 

9# 

10# Authors: 

11# - Colin Drieu 

12# - Dominique Chapelle 

13# - François Kimmig 

14# - Philippe Moireau 

15# 

16# PhysioBlocks is free software: you can redistribute it and/or modify it under the 

17# terms of the GNU Lesser General Public License as published by the Free Software 

18# Foundation, version 3 of the License. 

19# 

20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY 

21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 

22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along with 

25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>. 

26 

27from dataclasses import dataclass 

28from typing import Any 

29 

30import numpy as np 

31from numpy.typing import NDArray 

32 

33from physioblocks.computing import ( 

34 Expression, 

35 ModelComponent, 

36 Quantity, 

37 diff, 

38 mid_point, 

39) 

40from physioblocks.registers import register_type 

41from physioblocks.simulation import Time 

42 

43# Constant for the active law type id 

44ACTIVE_LAW_MACRO_HUXLEY_TWO_MOMENTS_TYPE_ID = "active_law_macro_huxley_two_moments" 

45 

46 

47@register_type(ACTIVE_LAW_MACRO_HUXLEY_TWO_MOMENTS_TYPE_ID) 

48@dataclass 

49class ActiveLawMacroscopicHuxleyTwoMoment(ModelComponent): 

50 r""" 

51 Describes the Macroscopic Huxley Two Moment implementation of the ative law. 

52 

53 **Internal Equations:** 

54 

55 .. math:: 

56 

57 \dot{K}_c + \bigl(|u|+\alpha|\dot{e}_c|\bigr) \, K_c - n_0(e_c) K_0 |u|_+ = 0 

58 

59 .. math:: 

60 

61 \dot{\lambda}_c + \frac{1}{2} \bigl(|u|+\alpha|\dot{e}_c|\bigr) \, \lambda_c - \sqrt{K_c} \, \dot{e}_c - \frac{n_0(e_c)}{\sqrt{K_c}} \biggl( T_0 - \frac{K_0\,\lambda_c}{2\sqrt{K_c}} \biggr) |u|_+ = 0 

62 

63 

64 .. math:: 

65 

66 \lambda_c - T_c/\sqrt{K_c} = 0 

67 

68 **Discretised:** 

69 

70 .. math:: 

71 

72 \frac{K_c^{n+1}-K_c^n}{\Delta t^n} + \Bigl(|u^{n+1}|+\alpha\Bigl|\frac{e_c^{n+1}-e_c^n}{\Delta t^n}\Bigr|\Bigr) \, K_c^{n+1} - n_0(e_c^n) K_0 \, |u^{n+1}|_+ = 0 

73 

74 .. math:: 

75 

76 \frac{\lambda_c^{n+1}-\lambda_c^n}{\Delta t^n} + \frac{1}{2}\Bigl(|u^{n+1}|+\alpha\Bigl|\frac{e_c^{n+1}-e_c^n}{\Delta t^n}\Bigr|\Bigr) \, \lambda_c^{n + \frac{1}{2}} - \sqrt{K_c^{n+1}} \, \frac{e_c^{n+1}-e_c^n}{\Delta t^n} - \frac{n_0(e_c^n)}{\sqrt{K_c^{n+1}}} \biggl( T_0 - \frac{K_0\,\lambda_c^{n + \frac{1}{2}}}{2\sqrt{K_c^{n+1}}} \biggr) |u^{n+1}|_+ = 0 

77 

78 .. math:: 

79 

80 T_c^{n + \frac{1}{2}\sharp} - \lambda_c^{n + \frac{1}{2}} \sqrt{K_c^{n+1}} = 0 

81 """ # noqa: E501 

82 

83 fib_deform: Quantity[np.float64] 

84 """:math:`e_c` the fiber deformation """ 

85 

86 active_stiffness: Quantity[np.float64] 

87 """:math:`K_c` the active stiffness""" 

88 

89 active_energy_sqrt: Quantity[np.float64] 

90 """:math:`\\lambda_c` the active energy sqrt""" 

91 

92 active_tension_discr: Quantity[np.float64] 

93 """:math:`T_c^{n + \\frac{1}{2}\\sharp}` the active tension discretisation""" 

94 

95 starling_abscissas: Quantity[NDArray[np.float64]] 

96 """Abscissas of :math:`n_0` the discretized starling function""" 

97 

98 starling_ordinates: Quantity[NDArray[np.float64]] 

99 """Ordinates of :math:`n_0` the discretize starling function""" 

100 

101 activation: Quantity[np.float64] 

102 """:math:`u` the activation function""" 

103 

104 crossbridge_stiffness: Quantity[np.float64] 

105 """:math:`K_0` the crossbridge stiffness""" 

106 

107 contractility: Quantity[np.float64] 

108 """:math:`T_0` the contractility""" 

109 

110 destruction_rate: Quantity[np.float64] 

111 """the destruction rate""" 

112 

113 time: Time 

114 """:math:`t` the simulation time""" 

115 

116 def active_law_residual(self) -> Any: 

117 """ 

118 Compute the residual of the active law. 

119 

120 :return: the residual value 

121 :rtype: np.array 

122 """ 

123 starling = np.interp( 

124 self.fib_deform.current, 

125 self.starling_abscissas.current, 

126 self.starling_ordinates.current, 

127 ) 

128 abs_plus_starling = np.clip(starling, 0.0, None) 

129 abs_plus_activation = np.clip(self.activation.new, 0.0, None) 

130 active_stiffness_dynamic_eq = ( 

131 (self.active_stiffness.new - self.active_stiffness.current) 

132 * self.time.inv_dt 

133 + self._get_deform_rate() * self.active_stiffness.new 

134 - abs_plus_starling 

135 * self.crossbridge_stiffness.current 

136 * abs_plus_activation 

137 ) 

138 active_energy_sqrt_dynamic_eq = ( 

139 diff(self.active_energy_sqrt) * self.time.inv_dt 

140 + 0.5 * self._get_deform_rate() * mid_point(self.active_energy_sqrt) 

141 - (abs_plus_starling / np.sqrt(self.active_stiffness.new)) 

142 * ( 

143 self.contractility.current 

144 - ( 

145 0.5 

146 * self.crossbridge_stiffness.current 

147 * mid_point(self.active_energy_sqrt) 

148 / (np.sqrt(self.active_stiffness.new)) 

149 ) 

150 ) 

151 * abs_plus_activation 

152 - np.sqrt(self.active_stiffness.new) 

153 * diff(self.fib_deform) 

154 * self.time.inv_dt 

155 ) 

156 

157 active_tension_discr_relation = self.active_tension_discr.new - np.sqrt( 

158 self.active_stiffness.new 

159 ) * mid_point(self.active_energy_sqrt) 

160 

161 return np.array( 

162 [ 

163 active_stiffness_dynamic_eq, 

164 active_energy_sqrt_dynamic_eq, 

165 active_tension_discr_relation, 

166 ], 

167 ) 

168 

169 def active_law_residual_dactive_stiffness(self) -> Any: 

170 """ 

171 Compute the residual partial derivative for the active stiffness 

172 

173 :return: the partial derivative 

174 :rtype: np.array 

175 """ 

176 starling = np.interp( 

177 self.fib_deform.current, 

178 self.starling_abscissas.current, 

179 self.starling_ordinates.current, 

180 ) 

181 abs_plus_starling = np.clip(starling, 0.0, None) 

182 abs_plus_activation = np.clip(self.activation.new, 0.0, None) 

183 

184 active_stiffness_active_law_dactive_stiffness = ( 

185 self.time.inv_dt + self._get_deform_rate() 

186 ) 

187 

188 active_energy_sqrt_active_law_dactive_stiffness = ( 

189 -abs_plus_activation 

190 * abs_plus_starling 

191 * ( 

192 -0.5 

193 * self.contractility.current 

194 * np.pow(self.active_stiffness.new, -1.5) 

195 + ( 

196 0.5 

197 * self.crossbridge_stiffness.current 

198 * mid_point(self.active_energy_sqrt) 

199 * np.pow(self.active_stiffness.new, -2) 

200 ) 

201 ) 

202 - 0.5 

203 * diff(self.fib_deform) 

204 * self.time.inv_dt 

205 / np.sqrt(self.active_stiffness.new) 

206 ) 

207 

208 active_tension_discr_relation_dactive_stiffness = ( 

209 -0.5 

210 * mid_point(self.active_energy_sqrt) 

211 / np.sqrt(self.active_stiffness.new) 

212 ) 

213 

214 return np.array( 

215 [ 

216 active_stiffness_active_law_dactive_stiffness, 

217 active_energy_sqrt_active_law_dactive_stiffness, 

218 active_tension_discr_relation_dactive_stiffness, 

219 ], 

220 ) 

221 

222 def active_law_residual_dactive_energy_sqrt(self) -> Any: 

223 """ 

224 Compute the residual partial derivative for the active energy sqrt 

225 

226 :return: the partial derivative 

227 :rtype: np.array 

228 """ 

229 starling = np.interp( 

230 self.fib_deform.current, 

231 self.starling_abscissas.current, 

232 self.starling_ordinates.current, 

233 ) 

234 abs_plus_starling = np.clip(starling, 0.0, None) 

235 abs_plus_activation = np.clip(self.activation.new, 0.0, None) 

236 

237 active_stiffness_active_law_dactive_energy_sqrt = 0.0 

238 

239 active_energy_sqrt_active_law_dactive_energy_sqrt = ( 

240 self.time.inv_dt 

241 + 0.25 

242 * ( 

243 abs(self.activation.new) 

244 + self.destruction_rate.current 

245 * (self.fib_deform.new - self.fib_deform.current) 

246 * self.time.inv_dt 

247 ) 

248 + 0.25 

249 * abs_plus_activation 

250 * abs_plus_starling 

251 * self.crossbridge_stiffness.current 

252 / self.active_stiffness.new 

253 ) 

254 

255 active_tension_discr_relation_dactive_energy_sqrt = -0.5 * np.sqrt( 

256 self.active_stiffness.new 

257 ) 

258 

259 return np.array( 

260 [ 

261 active_stiffness_active_law_dactive_energy_sqrt, 

262 active_energy_sqrt_active_law_dactive_energy_sqrt, 

263 active_tension_discr_relation_dactive_energy_sqrt, 

264 ], 

265 ) 

266 

267 def active_law_residual_dfib_deform(self) -> Any: 

268 """ 

269 Compute the residual partial derivative for the fiber deformation 

270 

271 :return: the partial derivative 

272 :rtype: np.array 

273 """ 

274 active_stiffness_active_law_dfib_deform = ( 

275 self.destruction_rate.current * self.active_stiffness.new * self.time.inv_dt 

276 ) 

277 

278 active_energy_sqrt_active_law_dfib_deform = -( 

279 -0.5 

280 * self.destruction_rate.current 

281 * mid_point(self.active_energy_sqrt) 

282 * self.time.inv_dt 

283 + np.sqrt(self.active_stiffness.new) * self.time.inv_dt 

284 ) 

285 

286 active_tension_discr_dfib_deform = 0.0 

287 

288 return np.array( 

289 [ 

290 active_stiffness_active_law_dfib_deform, 

291 active_energy_sqrt_active_law_dfib_deform, 

292 active_tension_discr_dfib_deform, 

293 ], 

294 ) 

295 

296 def active_law_residual_dactive_tension(self) -> Any: 

297 """ 

298 Compute the residual partial derivative for the active tension 

299 

300 :return: the partial derivative 

301 :rtype: np.array 

302 """ 

303 return np.array( 

304 [0.0, 0.0, 1.0], 

305 ) 

306 

307 def _get_deform_rate(self) -> Any: 

308 return np.fabs(self.activation.new) + ( 

309 self.destruction_rate.current 

310 * np.fabs(diff(self.fib_deform)) 

311 * self.time.inv_dt 

312 ) 

313 

314 

315# Define the expression for the residual of the ActiveLawMacroscopicHuxleyTwoMoment and 

316# its partial derivatives. 

317_active_law_macroscopic_huxley_two_moment_residual_expression = Expression( 

318 3, 

319 ActiveLawMacroscopicHuxleyTwoMoment.active_law_residual, 

320 { 

321 "active_stiffness": ActiveLawMacroscopicHuxleyTwoMoment.active_law_residual_dactive_stiffness, # noqa: E501 

322 "active_energy_sqrt": ActiveLawMacroscopicHuxleyTwoMoment.active_law_residual_dactive_energy_sqrt, # noqa: E501 

323 "active_tension_discr": ActiveLawMacroscopicHuxleyTwoMoment.active_law_residual_dactive_tension, # noqa: E501 

324 "fib_deform": ActiveLawMacroscopicHuxleyTwoMoment.active_law_residual_dfib_deform, # noqa: E501 

325 }, 

326) 

327 

328ActiveLawMacroscopicHuxleyTwoMoment.declares_internal_expression( 

329 "active_stiffness", 

330 _active_law_macroscopic_huxley_two_moment_residual_expression, 

331 1, 

332 0, 

333) 

334ActiveLawMacroscopicHuxleyTwoMoment.declares_internal_expression( 

335 "active_energy_sqrt", 

336 _active_law_macroscopic_huxley_two_moment_residual_expression, 

337 1, 

338 1, 

339) 

340ActiveLawMacroscopicHuxleyTwoMoment.declares_internal_expression( 

341 "active_tension_discr", 

342 _active_law_macroscopic_huxley_two_moment_residual_expression, 

343 1, 

344 2, 

345)