Skip to content

Commit

Permalink
Q10 implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@Ge0rges
Ge0rges committed Nov 13, 2024
1 parent 3a9771f commit 6068442
Show file tree
Hide file tree
Showing 8 changed files with 234 additions and 171 deletions.
Binary file removed Results/CB1 - Min \mu_{max} - Low m - Calculated \gamma_{cell}_sa.tif
View file
Binary file not shown.
Binary file removed Results/SA_result_object
View file
Binary file not shown.
30 changes: 6 additions & 24 deletions Results/values.csv
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,7 @@
Analysis title,Surrounding pOC,Brine pOC,Predicted brine pOC,Surrounding dOC,Brine dOC,Predicted brine dOC,Present cell density,Predicted end cell density,EEA Lower,EEA Upper,Real EEA,EEA predicted timespan,dOC/cell,Maintenance energy lower bounds,Maintenance energy upper bound,Simulation growth rate,Minimum growth rate,Minimum doubling time,Growth yield,Brine expansion factor
CB1 - Min $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.64e+13,1.49e+11,6.25e+12,4.08e+10,1.23e+12,6.5e+12,5.7e+06,5.71e+06,1.95e-01,8.02e-01,4.99e-01,-2.74e-03,1.57e+01,1.8071e-01,7.4366e-01,2.7673e-07,2.7673e-07,6.86e+03,2.4042e-05,3.33e-01
CB4 - Min $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.69e+13,4.97e+10,6.67e+12,7.4e+11,1.02e+12,7.48e+12,1.14e+07,1.14e+07,1.01e-01,4.83e-01,2.92e-01,-2.74e-03,5.40e+01,9.9397e-02,4.7491e-01,3.2418e-07,3.2418e-07,5.85e+03,1.7622e-04,3.33e-01
CBIW - Min $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.64e+13,2.38e+10,7.55e+12,4.08e+10,3.60e+11,6.77e+12,1.39e+08,1.24e+08,8.06e-03,5.84e-02,3.32e-02,-2.74e-03,1.57e+01,7.9006e-03,5.7210e-02,4.9535e-07,4.9535e-07,3.83e+03,9.8339e-04,3.33e-01
CB1 - Min $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,1.22e-02,8.12e+04,1.57e+01,1.8071e-01,7.4366e-01,2.7673e-07,2.7673e-07,6.86e+03,2.4042e-05,3.33e-01
CB4 - Min $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.68e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,1.22e-02,7.10e+04,5.40e+01,9.9397e-02,4.7491e-01,3.2418e-07,3.2418e-07,5.85e+03,1.7622e-04,3.33e-01
CBIW - Min $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,1.31e+13,4.08e+10,3.60e+11,1.18e+12,1.39e+08,1.24e+08,8.06e-03,5.84e-02,1.22e-02,4.87e+04,1.57e+01,7.9006e-03,5.7210e-02,4.9535e-07,4.9535e-07,3.83e+03,9.8339e-04,3.33e-01
CB1 - Min $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.64e+13,1.49e+11,1.63e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,4.99e-01,-2.74e-03,1.57e+01,1.8071e-01,7.4366e-01,2.7673e-07,2.7673e-07,6.86e+03,5.8423e-06,3.33e-01
CB4 - Min $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.69e+13,4.97e+10,1.57e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,2.92e-01,-2.74e-03,5.40e+01,9.9397e-02,4.7491e-01,3.2418e-07,3.2418e-07,5.85e+03,3.6886e-05,3.33e-01
CBIW - Min $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.64e+13,2.38e+10,1.63e+13,4.08e+10,3.60e+11,7.02e+08,1.39e+08,7.35e+00,8.06e-03,5.84e-02,3.32e-02,-2.74e-03,1.57e+01,7.9006e-03,5.7210e-02,4.9535e-07,4.9535e-07,3.83e+03,1.3592e-04,3.33e-01
CB1 - Min $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,1.22e-02,8.12e+04,1.57e+01,1.8071e-01,7.4366e-01,2.7673e-07,2.7673e-07,6.86e+03,5.8423e-06,3.33e-01
CB4 - Min $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.69e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,1.22e-02,7.10e+04,5.40e+01,9.9397e-02,4.7491e-01,3.2418e-07,3.2418e-07,5.85e+03,3.6886e-05,3.33e-01
CBIW - Min $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,1.64e+13,4.08e+10,3.60e+11,7.02e+08,1.39e+08,7.35e+00,8.06e-03,5.84e-02,1.22e-02,4.87e+04,1.57e+01,7.9006e-03,5.7210e-02,4.9535e-07,4.9535e-07,3.83e+03,1.3592e-04,3.33e-01
CB1 - Max $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.64e+13,1.49e+11,0.e+00,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,4.99e-01,-2.74e-03,1.57e+01,1.8071e-01,7.4366e-01,6.e-02,2.7673e-07,6.86e+03,8.3904e-01,3.33e-01
CB4 - Max $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.69e+13,4.97e+10,0.e+00,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,2.92e-01,-2.74e-03,5.40e+01,9.9397e-02,4.7491e-01,1.6e-02,3.2418e-07,5.85e+03,8.969e-01,3.33e-01
CBIW - Max $\mu_{max}$ - Low $m$ - Calculated $\gamma_{cell}$,1.64e+13,2.38e+10,0.e+00,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,5.84e-02,3.32e-02,-2.74e-03,1.57e+01,7.9006e-03,5.7210e-02,6.e-02,4.9535e-07,3.83e+03,9.9168e-01,3.33e-01
CB1 - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,1.22e-02,8.12e+04,1.57e+01,1.8071e-01,7.4366e-01,6.e-02,2.7673e-07,6.86e+03,8.3904e-01,3.33e-01
CB4 - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.68e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,1.22e-02,7.10e+04,5.40e+01,9.9397e-02,4.7491e-01,1.6e-02,3.2418e-07,5.85e+03,8.969e-01,3.33e-01
CBIW - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,0.e+00,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,5.84e-02,1.22e-02,4.87e+04,1.57e+01,7.9006e-03,5.7210e-02,6.e-02,4.9535e-07,3.83e+03,9.9168e-01,3.33e-01
CB1 - Max $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.64e+13,1.49e+11,1.63e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,4.99e-01,-2.74e-03,1.57e+01,1.8071e-01,7.4366e-01,6.e-02,2.7673e-07,6.86e+03,5.5883e-01,3.33e-01
CB4 - Max $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.69e+13,4.97e+10,1.57e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,2.92e-01,-2.74e-03,5.40e+01,9.9397e-02,4.7491e-01,1.6e-02,3.2418e-07,5.85e+03,6.4547e-01,3.33e-01
CBIW - Max $\mu_{max}$ - High $m$ - Calculated $\gamma_{cell}$,1.64e+13,2.38e+10,1.63e+13,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,5.84e-02,3.32e-02,-2.74e-03,1.57e+01,7.9006e-03,5.7210e-02,6.e-02,4.9535e-07,3.83e+03,9.4274e-01,3.33e-01
CB1 - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,8.02e-01,1.22e-02,8.12e+04,1.57e+01,1.8071e-01,7.4366e-01,6.e-02,2.7673e-07,6.86e+03,5.5883e-01,3.33e-01
CB4 - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.69e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,4.83e-01,1.22e-02,7.10e+04,5.40e+01,9.9397e-02,4.7491e-01,1.6e-02,3.2418e-07,5.85e+03,6.4547e-01,3.33e-01
CBIW - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,1.64e+13,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,5.84e-02,1.22e-02,4.87e+04,1.57e+01,7.9006e-03,5.7210e-02,6.e-02,4.9535e-07,3.83e+03,9.4274e-01,3.33e-01
CB1 - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,3.03e+00,1.22e-02,5.42e+04,1.57e+01,1.8071e-01,2.8111e+00,6.e-02,1.0647e-06,1.78e+03,8.3904e-01,3.33e-01
CB4 - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.68e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,1.64e+00,1.22e-02,5.21e+04,5.40e+01,9.9397e-02,1.6151e+00,1.6e-02,1.1122e-06,1.71e+03,8.969e-01,3.33e-01
CBIW - Max $\mu_{max}$ - Low $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,0.e+00,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,1.51e-01,1.22e-02,4.54e+04,1.57e+01,7.9006e-03,1.4812e-01,6.e-02,1.2834e-06,1.48e+03,9.9168e-01,3.33e-01
CB1 - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,1.49e+11,1.64e+13,4.08e+10,1.23e+12,0.e+00,5.7e+06,0.e+00,1.95e-01,3.03e+00,1.22e-02,5.42e+04,1.57e+01,1.8071e-01,2.8111e+00,6.e-02,1.0647e-06,1.78e+03,2.5099e-01,3.33e-01
CB4 - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.69e+13,4.97e+10,1.69e+13,7.4e+11,1.02e+12,0.e+00,1.14e+07,0.e+00,1.01e-01,1.64e+00,1.22e-02,5.21e+04,5.40e+01,9.9397e-02,1.6151e+00,1.6e-02,1.1122e-06,1.71e+03,3.4869e-01,3.33e-01
CBIW - Max $\mu_{max}$ - High $m$ - Measured $\gamma_{cell}$,1.64e+13,2.38e+10,1.64e+13,4.08e+10,3.60e+11,0.e+00,1.39e+08,0.e+00,8.06e-03,1.51e-01,1.22e-02,4.54e+04,1.57e+01,7.9006e-03,1.4812e-01,6.e-02,1.2834e-06,1.48e+03,8.6413e-01,3.33e-01
3 changes: 1 addition & 2 deletions analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,8 +119,7 @@ def run_analysis(self, do_sensitivity_analysis: bool, cached_SA: bool = False):
print("do_SA and cached_SA were both True. Loading results from cache only.")

elif do_sensitivity_analysis:
self.sensitivity_analysis_result = run_sensitivity_analysis(
self.scenario) if do_sensitivity_analysis else None
self.sensitivity_analysis_result = run_sensitivity_analysis(self.scenario) if do_sensitivity_analysis else None


def estimate_me_bounds(scenario: Scenario):
Expand Down
47 changes: 24 additions & 23 deletions main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,38 +93,39 @@ def log_results(analyses, cached_SA):

# Plot for all scenarios all analyses
all_analyses_fig = plot_all_scenarios_all_analyses(analyses)
single_datatype_plot_cells = plot_one_result_type_all_analyses(analyses, "Cells", 0)
single_datatype_plot_cells = plot_one_result_type_all_analyses(analyses, "Cells")

all_analyses_fig.savefig("Results/all_model_outputs.tif", format="tif", dpi=300)
single_datatype_plot_cells.savefig("Results/all_model_outputs_just_cells.tif", format="tif", dpi=300)

# Example usage
interactive_plot = plot_one_result_type_all_analyses_interactive_plotly(all_analyses, "Cells")
interactive_plot.show()

return


if __name__ == "__main__": # Generates all figures and data points.
# All sensitivity analyses should be the same.
scenarios = [cb1_scenario(), cb4_scenario(), cbiw_scenario()]
all_analyses = []

# On every scenario, try every analysis configuration. Run a sensitivity analysis only on the first one.
do_SA = False # If true, run SA. If do_SA_from_file is also true, results will be taken from cached file.
cached_SA = not do_SA # If true, loads SA results from file and then plots it.
for use_minimum_growth_rate in [True, False]:
for use_me_lower_bound in [True, False]:
for use_eea_average in [True, False]:
for scenario in scenarios:
a = Analysis(scenario, use_minimum_growth_rate, use_me_lower_bound, use_eea_average)
cached_results = False
if not cached_results:
# All sensitivity analyses should be the same.
scenarios = [cb1_scenario(), cb4_scenario(), cbiw_scenario()]
all_analyses = []

# Run SA on first config pair - first scenario, only.
a.run_analysis(do_sensitivity_analysis=do_SA, cached_SA=cached_SA)
do_SA = False
cached_SA = False
# On every scenario, try every analysis configuration.
for use_me_lower_bound in [True, False]:
for scenario in scenarios:
a = Analysis(scenario, use_minimum_growth_rate=False, use_me_lower_bound=use_me_lower_bound, use_eea_average=False)
a.run_analysis(do_sensitivity_analysis=False, cached_SA=False)
all_analyses.append(a)

all_analyses.append(a)
# Save plots and values of all results.
with open("Results/all_analyses", "wb") as aa_results_file:
pickle.dump(all_analyses, aa_results_file)

# Save plots and values of all results.
log_results(all_analyses, cached_SA)
log_results(all_analyses, False)

# Generate the hypothetical growth scenarios figure.
hg_fig = hypothetical_growth_scenarios()
hg_fig.savefig("Results/" + "hypothetical_growth.tif", format="tif", dpi=300)
else: # Load results and create plots
with open("Results/all_analyses", "rb") as aa_results_file:
all_analyses = pickle.load(aa_results_file)
log_results(all_analyses, False)
40 changes: 35 additions & 5 deletions model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,20 @@ using ForwardDiff
using GlobalSensitivity
using Statistics
using Printf
using Interpolations
using CSV
using DataFrames


# Create the linear interpolation object for temperature catch. Do it here once.
time_df = sort!(CSV.read("temp_reconstruction.csv", DataFrame)) # Load the CSV file into a DataFrame. CSV is time, temperature
temp_interp = LinearInterpolation(time_df[:,1], time_df[:,2]) # Create an interpolation object

# Some variables
Q10 = 1.5 # Assumed Q10 value for EEA and growth rate
mu_T0 = -10 # Original temperature of the max growth rate.
eea_T0 = -6 # Original temperature of the eea rate.
println("Q10 adjustments in effect. Constants defined in top of model.jl file.")


# Runs the model using solve_model and packages results nicely.
Expand All @@ -15,7 +29,7 @@ function run_model(p, u0)
u0 = convert(Array{Float64}, u0)

prob = ODEProblem(model, u0[1:end-1], (0.0, last(u0)), p)
sol = solve(prob, Rosenbrock23(), callback=cbs, maxiters=1e6)
sol = solve(prob, Rosenbrock23(), callback=cbs, maxiters=1e7, tstops=time_df[:,1])

# Return the solution array - [pOC, dOC, DIC, Cells, t]
return [[x[1] for x in sol.u], [x[2] for x in sol.u], [x[3] for x in sol.u], [x[4] for x in sol.u], sol.t]
Expand Down Expand Up @@ -84,8 +98,12 @@ function make_callbacks(;sensitivity_analysis=false, p=nothing)
do_nothing(integrator) = nothing
max_cb = ContinuousCallback(max_condition, do_nothing)

# Growth rate condition
growth_condition(u, t, integer) = (temp_interp(t) < -15)
gr_cb = DiscreteCallback(growth_condition, do_nothing)

# Min condition for EEA rate
eea_min_condition(u, t, integrator) = integrator.p[10] * u[4] - u[1]
eea_min_condition(u, t, integrator) = integrator.p[10] * Q10^((temp_interp(t) - eea_T0) / 10) * u[4] - u[1]
eea_min_cb = ContinuousCallback(eea_min_condition, do_nothing)

# Min condition for DIC fixation rate
Expand All @@ -104,7 +122,7 @@ function make_callbacks(;sensitivity_analysis=false, p=nothing)
zero_cb = DiscreteCallback(zero_condition, zero_out!)

# Callback list
return CallbackSet(carbon_add_cb, max_cb, eea_min_cb, ic_min_cb, zero_cb, PositiveDomain())
return CallbackSet(carbon_add_cb, max_cb, eea_min_cb, gr_cb, ic_min_cb, zero_cb, PositiveDomain())
end


Expand All @@ -131,8 +149,18 @@ function model(du, u, p, t)
cell_count = u[4]

## CELL COUNT

# Adjust growth rate for Q10
current_temp = temp_interp(t) # Find temperature at current time
if current_temp < -15
# Set growth rate to 0 if temperature is below -15°C
adjusted_growth_rate = 0
else
adjusted_growth_rate = mu_max * Q10^((current_temp - mu_T0) / 10)
end

# Growth
growth = mu_max * (dOC_content / (Kd + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)
growth = adjusted_growth_rate * (dOC_content / (Kd + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

# Organic carbon requirement
required_dOC_per_cell = maintenance_per_cell
Expand All @@ -150,7 +178,9 @@ function model(du, u, p, t)
fixed_carbon = min(inorganic_carbon_fixing_rate, inorganic_carbon_content)

# EEA rate
eea_removal = min(eea_rate*cell_count, pOC_content)
adjusted_eea_rate = eea_rate * Q10^((current_temp - eea_T0) / 10)

eea_removal = min(adjusted_eea_rate*cell_count, pOC_content)

# Particulate Organic carbon
du[1] = pOC_input_rate - eea_removal
Expand Down
Loading

0 comments on commit 6068442

Please sign in to comment.