StemCellTheoryPaperCode/fig_3abc_plot_two_compartment_intro.py at master · jvzonlab/StemCellTheoryPaperCode · GitHub

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
"""Used to give a general overview of the two compartment model"""
import math
from collections import OrderedDict
from typing import Dict, List, Tuple

import matplotlib
import numpy
from matplotlib.axes import Axes
from matplotlib.colorbar import Colorbar

from stem_cell_model import sweeper, tools, two_compartment_model
import matplotlib.pyplot as plt

from stem_cell_model.parameters import SimulationConfig, SimulationParameters
from stem_cell_model.results import SimulationResults

matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['svg.fonttype'] = 'none'

STEPS_ALONG_AXIS = 40
COLOR_MAP = "gnuplot"


def plot_alpha_m_against_alpha_n_and_phi_at_1():
    """Builds an image."""
    image = numpy.full((STEPS_ALONG_AXIS, STEPS_ALONG_AXIS), numpy.nan, dtype=numpy.float64)
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.phi[0] == 1:
            image_x = int(-params.alpha[1] * (STEPS_ALONG_AXIS - 1))
            image_y = int(params.alpha[0] * (STEPS_ALONG_AXIS - 1))

            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            image[image_y, image_x] = statistics.d_coeff_var

    return image


def plot_phi_against_opposite_alpha():
    """Builds an image."""
    image = numpy.full((STEPS_ALONG_AXIS, STEPS_ALONG_AXIS), numpy.nan, dtype=numpy.float64)
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.alpha[0] == -params.alpha[1]:
            image_y = int(params.phi[0] * (STEPS_ALONG_AXIS - 1))
            image_x = int(params.alpha[0] * (STEPS_ALONG_AXIS - 1))
            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            image[image_y, image_x] = statistics.d_coeff_var

    return image


def get_example_line_for_phi_1_opposite_alpha() -> Tuple[List[float], List[float]]:
    """Line along the largest variation."""
    values = dict()
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.alpha[0] == -params.alpha[1] and params.phi[0] == 1:
            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            values[params.alpha[0]] = statistics.d_coeff_var

    x_values = list()
    y_values = list()
    for key, value in sorted(values.items()):
        x_values.append(key)
        y_values.append(value)
    return x_values, y_values


def get_example_line_for_opposite_alpha_and_matching_phi() -> Tuple[List[float], List[float]]:
    """Line along the lowest phi for every alpha. This is one of the optimal lines for minimizing fluctuations.."""
    values = dict()
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.alpha[0] == -params.alpha[1] and params.phi[0] == params.alpha[0]:
            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            values[params.alpha[0]] = statistics.d_coeff_var

    x_values = list()
    y_values = list()
    for key, value in sorted(values.items()):
        x_values.append(key)
        y_values.append(value)
    return x_values, y_values


def get_example_line_for_opposite_alpha_and_higher_phi() -> Tuple[List[float], List[float]]:
    """Line along (0.5 * lowest phi) + 0.5 for every alpha. So we're halfway the lowest phi
    and the maximal phi (which is always 1)."""
    values = dict()
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.alpha[0] == -params.alpha[1] and (0.5 * params.alpha[0]) + 0.5 == params.phi[0]:
            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            values[params.alpha[0]] = statistics.d_coeff_var

    x_values = list()
    y_values = list()
    for key, value in sorted(values.items()):
        x_values.append(key)
        y_values.append(value)
    return x_values, y_values


def get_second_example_line_for_phi_1() -> Tuple[List[float], List[float]]:
    """The line for phi=1 and a_n = -a_m + 0.5."""
    values = dict()
    for params, multi_run_stats in sweeper.load_sweep_results("two_comp_sweep_data_fixed_D"):
        if params.alpha[0] == -params.alpha[1] + 0.5 and params.phi[0] == 1:
            statistics = tools.get_single_parameter_set_statistics(multi_run_stats)
            values[params.alpha[0]] = statistics.d_coeff_var

    x_values = list()
    y_values = list()
    for key, value in sorted(values.items()):
        x_values.append(key)
        y_values.append(value)
    return x_values, y_values


fig, (ax_a, ax_b, ax_c) = plt.subplots(nrows=3,figsize=(4.181, 6.498),
                  gridspec_kw={"height_ratios": [0.47, 0.47, 0.06]})

image_phi1 = plot_alpha_m_against_alpha_n_and_phi_at_1()
image_varphi = plot_phi_against_opposite_alpha()
images_max = max(numpy.nanmax(image_varphi), numpy.nanmax(image_phi1))  # To keep the color scale of both images the same
images_max = math.ceil(images_max)


# Draw the image for phi =1
ax_a.set_title("Coeff of var $D(t)$ at $\phi = 1$")
ax_a.set_facecolor("#b2bec3")
ax_a.imshow(image_phi1, extent=(0, -1, 1, 0), aspect=1, cmap=COLOR_MAP, interpolation="nearest", vmin=0, vmax=images_max)
ax_a.set_xlabel("$\\alpha_m$")
ax_a.set_ylabel("$\\alpha_n$")

# Draw the image
ax_b.set_title("Coeff of var $D(t)$")
ax_b.set_facecolor("#b2bec3")
ax_b_image = ax_b.imshow(image_varphi, extent=(0, 1, 1, 0), aspect=1, cmap=COLOR_MAP, interpolation="nearest", vmin=0, vmax=images_max)
ax_b.invert_yaxis()
ax_b.set_xlabel("$\\alpha_n, -\\alpha_m$")
ax_b.set_ylabel("$\phi$")

colorbar: Colorbar = fig.colorbar(ax_b_image, cax=ax_c, orientation="horizontal")
fig.tight_layout()
plt.show()

# Now draw values along lines
fig, (ax_a, ax_b) = plt.subplots(nrows=2,figsize=(4.181, 4.498), sharex="col", sharey="col")
ax_a.plot(*get_example_line_for_phi_1_opposite_alpha(), color="#81ecec", label="$\\phi=1,\\alpha_n=-\\alpha_m$")
ax_a.plot(*get_second_example_line_for_phi_1(), color="#fd79a8", label="$\\phi=1,\\alpha_n=-\\alpha_m + 0.5$")
ax_a.set_ylabel("Coeff var D")
ax_a.set_xlabel("$\\alpha_n$")
ax_a.legend()

ax_b.plot(*get_example_line_for_opposite_alpha_and_matching_phi(), color="#d63031", label="$\\phi=\\alpha_n$")
ax_b.plot(*get_example_line_for_opposite_alpha_and_higher_phi(), color="#000000", label="$\\phi=0.5*\\alpha_n + 0.5$")
ax_b.set_ylabel("Coeff var D")
ax_b.set_xlabel("$\\alpha_n, -\\alpha_m$")
ax_b.legend()

ax_b.set_xlim(0, 1)
ax_b.set_ylim(0, 1)

plt.tight_layout()
plt.show()


# Now draw examples
_, (ax_left, ax_left_histogram, ax_right, ax_right_histogram) = plt.subplots(nrows=1, ncols=4, sharey="all",
                                                                             gridspec_kw={'width_ratios': [3, 1, 3, 1]})
random = numpy.random.Generator(numpy.random.MT19937(seed=1))
config = SimulationConfig(t_sim=750, random=random, track_n_vs_t=True)
T = (16.153070175438597, 3.2357834505600382)  # Based on measured values


def _plot_line(ax: Axes, results: SimulationResults):
    ax.plot(results.n_vs_t[:, 0], results.n_vs_t[:, 1] + results.n_vs_t[:, 2], color="black", alpha=0.5,
                 linewidth=2)
    if results.n_vs_t[-1, 1] + results.n_vs_t[-1, 2] == 0:  # Died
        ax.plot(results.n_vs_t[-1, 0], 0, "X", color="red")

def _plot_histogram(ax: Axes, config: SimulationConfig, params: SimulationParameters):
    bins = numpy.arange(0, 60, 1)
    counts = list()

    for i in range(2000):
        if i % 100 == 0:
            print(i)
        results = two_compartment_model.run_simulation(config, params)
        last_count = results.n_vs_t[-1, 1] + results.n_vs_t[-1, 2]
        counts.append(last_count)
    ax.hist(counts, orientation="horizontal", bins=bins, color="black")
    ax.axis("off")

# Left: low phi, low alpha
params = SimulationParameters.for_D_alpha_and_phi(D=30, phi=0.05, T=T, alpha_n=0, alpha_m=0)
for _ in range(6):
    results = two_compartment_model.run_simulation(config, params)
    _plot_line(ax_left, results)
_plot_histogram(ax_left_histogram, config, params)

# Right: large phi, large alpha
params = SimulationParameters.for_D_alpha_and_phi(D=30, phi=0.95, T=T, alpha_n=0.95, alpha_m=-0.95)
for _ in range(6):
    results = two_compartment_model.run_simulation(config, params)
    _plot_line(ax_right, results)
_plot_histogram(ax_right_histogram, config, params)

ax_left.set_ylabel("Proliferating cells")
ax_left.set_xlabel("Time (h)")
ax_left.set_ylim(0, 55)  # Propagates to other axis because of sharey

plt.show()