Commit dbe41816 authored by Félix Hartmann's avatar Félix Hartmann
Browse files

[feature] The coef of polar transport can be controled by local

concentration of signals.
parent bdefba2f
......@@ -63,7 +63,7 @@ class CellFile(object):
If True, strains are inversely proportional to cell diameters.
data_time_unit : int
time unit of data, in seconds, including signals (but not
diffusion constant)
diffusion constant and time steps)
date : datetime object
date of initialization of the file
"""
......@@ -651,9 +651,14 @@ class CellFile(object):
M_i[self.lower_crossed_mask] += dt * P / RdX[1:]
elif mechanism == "simple unidirectional active transport":
RdX_flux_i = 1/self.D_real[self.flux_cells]
p, q, d = transport.p, transport.q, transport.d
N = self.nb_cells
M_i = np.zeros((N, N))
p, q, d = transport.p, transport.q, transport.d
p = p * np.ones(N)
if transport.signals_controling_polar_transport:
# Signal concentrations are multiplied together.
conc = self.cell_conc[transport.polar_control_indices]
p += transport.polar_prefactor * np.prod(conc, axis=0)
if transport.only_in_cambial_cells:
polar_zone = self.zone_1.copy()
polar_zone_lower = polar_zone.copy()
......@@ -682,17 +687,19 @@ class CellFile(object):
Dupper = q * dt / self.D_real[:-1]
# additional terms for polar transport
D[polar_zone] -= \
p * dt / self.D_real[polar_zone]
p[polar_zone] * dt / self.D_real[polar_zone]
if transport.carrier_orientation == "toward the xylem":
Dlower[polar_zone_lower] += \
p * dt / self.D_real[polar_zone_lower + 1]
p[polar_zone_lower] * dt \
/ self.D_real[polar_zone_lower + 1]
if N-1 in polar_zone: # no polar transport at boundary
D[-1] += p * dt / self.D_real[-1]
if transport.carrier_orientation == "toward the cambium":
D[-1] += p[-1] * dt / self.D_real[-1]
elif transport.carrier_orientation == "toward the cambium":
Dupper[polar_zone_upper - 1] += \
p * dt / self.D_real[polar_zone_upper - 1]
p[polar_zone_upper] * dt \
/ self.D_real[polar_zone_upper - 1]
if 0 in polar_zone: # no polar transport at boundary
D[0] += p * dt / self.D_real[0]
D[0] += p[0] * dt / self.D_real[0]
np.fill_diagonal(M_i, D)
np.fill_diagonal(M_i[1:, :-1], Dlower)
np.fill_diagonal(M_i[:-1, 1:], Dupper)
......@@ -759,19 +766,21 @@ class CellFile(object):
"""
if activation_time is None:
activation_delay = start_time
activation_delay = start_time.total_seconds()
else:
activation_delay = (activation_time - start_time).total_seconds()
total_time = (stop_time - start_time).total_seconds()
self.start_time, self.stop_time = start_time, stop_time
nb_input_samples = \
int((total_time + snapshot_time_step) / input_time_step) + 2
self.date_list = self.date_list[:1]
self.growth_time_step = growth_time_step
self.conversion = self.data_time_unit / snapshot_time_step
self.initial_nodes_per_cell = nodes_per_cell
self.nodes_per_cell = nodes_per_cell*np.ones(self.nb_cells, dtype=int)
self.nb_nodes = self.nodes_per_cell.sum()
self.conc_operator = self.compute_concentration_operator(
self.nodes_per_cell)
self.conc_operator = \
self.compute_concentration_operator(self.nodes_per_cell)
# Construct the list of signals
signals = [division_identity_signal] + [enlargement_identity_signal] \
+ list(control_signals)
......@@ -783,6 +792,19 @@ class CellFile(object):
self.enl_id_index = self.signals.index(enlargement_identity_signal)
self.control_indices = [i for i, sig in enumerate(self.signals)
if sig in control_signals]
# Set values associated with the signals
for signal in self.signals:
signal.set_values(start_time.days,
start_time.days + (nb_input_samples - 1) *
input_time_step / self.data_time_unit,
nb_input_samples,
self.signals)
# Signal transport mechanisms
self.transport_mechanisms = [signal.transport_mechanism
for signal in self.signals]
self.transports = [signal.transport for signal in self.signals]
# List of time steps dt (one for each signal)
self.dt = np.zeros(self.nb_signals)
# length of the reference domaine Omega_0
self.L0 = self.nb_cells * self.initial_CRD
# mesh internode in each cell
......@@ -795,11 +817,6 @@ class CellFile(object):
# cumulative cell growth (R)
self.R = np.exp(self.int_S)
self.D_real = self.R * self.D
# signal characteristics
self.transport_mechanisms = [signal.transport_mechanism
for signal in self.signals]
self.transports = [signal.transport for signal in self.signals]
self.dt = np.zeros(self.nb_signals)
# empty concentration arrays
u, PIN = [], []
for mechanism in self.transport_mechanisms:
......@@ -814,8 +831,6 @@ class CellFile(object):
u.append(u_i)
PIN.append(PIN_i)
# boundary conditions
nb_input_samples = \
self.signals[0].boundary_conditions()[0].values.shape[0]
self.BC_values = [np.zeros((2, nb_input_samples))
for i in range(self.nb_signals)]
self.conc_indices = []
......@@ -851,12 +866,11 @@ class CellFile(object):
# initialize imposed concentrations
u[i][conc_indices_i] = self.BC_values[i][conc_indices_i][:, 0]
if transport_mechanism == "uniform diffusion":
if signal.transport.steady_state:
if transport.steady_state:
self.dt[i] = input_time_step
else:
# time step given by the Courant–Friedrichs–Lewy condition
self.dt[i] = self.dX_cell.min()**2 \
/ (2 * signal.transport.D)
self.dt[i] = self.dX_cell.min()**2 / (2 * transport.D)
# set concentration boundary conditions,
# with linear interpolation within the domain
if signal.nb_concentrations() == 2:
......@@ -866,7 +880,7 @@ class CellFile(object):
u[i] = np.linspace(u0, uL, self.nb_nodes)
elif transport_mechanism \
== "simple unidirectional active transport":
p, q, d = transport.p, transport.q, transport.d
p, q, d = transport.p_max(), transport.q, transport.d
p = abs(p)
Dmin = self.D_real.min()
self.dt[i] = 0.95 * Dmin / (p + 2*q + d*Dmin)
......@@ -875,17 +889,17 @@ class CellFile(object):
u[i][conc_indices_i] = conc
elif transport_mechanism == "Smith model (2006)":
self.dt[i] = 0.005
u[i], PIN[i] = signal.transport.evolution_conc(
u[i], PIN[i], self.D_real, self.S, (None, None, None),
self.dt[i], 5000)
u[i], PIN[i] = transport.evolution_conc(
u[i], PIN[i], self.D_real, self.S, (None, None, None),
self.dt[i], 5000)
elif transport_mechanism == "imposed cell concentrations":
phloem_BC, xylem_BC = None, None
if 0 in self.conc_indices[i]:
phloem_BC = self.BC_values[i][0][0]
if -1 in self.conc_indices[i]:
xylem_BC = self.BC_values[i][1][0]
u[i] = signal.transport.cell_concentrations(
0, self.nb_cells, phloem_BC, xylem_BC)
u[i] = transport.cell_concentrations(
0, self.nb_cells, phloem_BC, xylem_BC)
self.set_conc_and_flux_masks(self.conc_indices, self.flux_indices)
......
......@@ -168,17 +168,6 @@ def growth(simu):
start_time = timedelta(days=simu.start_time)
stop_time = timedelta(days=simu.stop_time)
activation_time = timedelta(days=simu.activation_time)
total_time = (stop_time - start_time).total_seconds()
nb_input_samples = int((total_time + simu.snapshot_time_step)
/ simu.input_time_step) + 2
# precomputation of time-dependent values
signals = simu.signals
for signal in signals:
signal.set_time_dependant_values(
simu.start_time,
simu.start_time + (nb_input_samples - 1) *
simu.input_time_step / time_units["day"],
nb_input_samples)
simu.cell_file = CellFile(
nb_cells=simu.nb_cells,
initial_CRD=simu.cell_diameter,
......@@ -542,8 +531,8 @@ class ControlPanel(HasTraits):
style='custom',
show_label=False),
HGroup(
Item('save_settings', show_label=False),
Item('load_settings', show_label=False)
Item('load_settings', show_label=False),
Item('save_settings', show_label=False)
),
label="Simulation",
dock='tab'
......
......@@ -87,15 +87,27 @@ class Morphogen(HasTraits):
nb_flux += 1
return nb_flux
def set_time_dependant_values(self, start, stop, num):
"""Set the values taken by the boundary conditions during the
simulation.
def set_values(self, start, stop, num, signals=None):
"""
Set various values, time-variable or not, associated with the signal.
Namely:
- Values taken by the boundary conditions during the simulation.
- Values taken by time-dependant parameters of the transport mechanism.
- Information on how signal transport depends on other signals.
Parameters
----------
start, stop and num: ints
Parameters for creating time points with Numpy's linspace function.
signals: iterable of Morphogen instances, optional
All signals involved in the simulation.
"""
# boundary conditions values
for bcond in self.boundary_conditions():
bcond.set_values(start, stop, num)
# values of the time-dependant parameters of the transport mechanism
self.transport.set_values(start, stop, num)
# and information on how signal transport depends on other signals
self.transport.set_values(start, stop, num, signals=signals)
def __repr__(self):
return self.name
......
......@@ -10,7 +10,7 @@ from __future__ import (unicode_literals, absolute_import,
import numpy as np
from traits.api import (HasTraits, Instance, Bool, Int, Float, Enum, Code)
from traits.api import (HasTraits, Instance, Bool, Int, Float, Str, Enum, Code)
from traitsui.api import (View, Item, Group, Handler, CodeEditor)
from evolution_function import (FunctionHandler, EvolutionFunction, Constant,
......@@ -47,7 +47,7 @@ class Transport(HasTraits):
"""
An empty class from which all particular transport type classes inherit.
"""
def set_values(self, start, stop, num):
def set_values(self, start, stop, num, signals=None):
pass
......@@ -128,7 +128,7 @@ class UniformDiffusion(Transport):
label="Steady state approximation",
desc="if the steady state approximation can be made.")
def set_values(self, start, stop, num):
def set_values(self, start, stop, num, signals=None):
"""
Set the values of the time-varying parameters.
"""
......@@ -209,6 +209,9 @@ class SimpleUnidirectionalActiveTransport(Transport):
"""
Compute the values taken by the polar transport coefficient during the
simulation.
- start, stop and num are the parameters for creating time points with
Numpy's linspace function.
"""
X = np.linspace(start, stop, num)
return self.p_function.f(X)
......@@ -229,6 +232,10 @@ class SimpleUnidirectionalActiveTransport(Transport):
raise AttributeError("Sorry, the decay rate can not assigned that way."
"Use the p_function trait instead.")
def p_max(self):
"""Return the maximal value of p."""
return np.max(self.p_values)
q = Float(1,
label="Permeability rate (µm.s¯¹)",
desc="permeability rate of the membranes, in µm.s¯¹.")
......@@ -251,13 +258,49 @@ class SimpleUnidirectionalActiveTransport(Transport):
only_in_cambial_cells = Bool(
label="Polar transport only in cambial cells",
desc="if polar transport is limited to cambial cells.")
desc="whether polar transport is limited to cambial cells.")
signals_controling_polar_transport = Str(
"",
label="Signals controling_polar_transport",
desc="which signals, if any, control the value of the polar transport "
"coefficient. Coma-separated list of signal names. Signals "
"listed twice are raised to the second power, and so on.")
polar_prefactor = Float(
0,
label="Polar transport prefactor",
desc="prefactor k_p multiplying signal concentrations involve in the "
"value of the polar transport coefficient.")
polar_control_indices = []
def set_values(self, start, stop, num):
def set_values(self, start, stop, num, signals=None):
"""
Set the values of the time-varying parameters.
Set values associated with the signal.
Parameters
----------
start, stop and num: ints
Parameters for creating time points with Numpy's linspace function.
signals: iterable of Morphogen instances, optional
All signals involved in the simulation.
"""
# First, set the values of constitutive polar transport.
self.set_p_values(start, stop, num)
# Second, look for the indices of signal controling the coefficient of
# polar transport.
if signals is not None:
names = [sig.strip() for sig in
self.signals_controling_polar_transport.split(",") if sig]
self.polar_control_indices = []
for name in names:
try:
signal_names = [signal.name for signal in signals]
self.polar_control_indices.append(signal_names.index(name))
except ValueError:
print("Warning: {0} is not a valid signal name for polar "
"transport control".format(name))
general_group = Group(Item('d'),
Item('q'),
......@@ -286,6 +329,8 @@ class SimpleUnidirectionalActiveTransport(Transport):
# Item('s'),
Item('carrier_orientation'),
Item('only_in_cambial_cells'),
Item('signals_controling_polar_transport'),
Item('polar_prefactor'),
springy=True
),
handler=FunctionHandler(),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment