### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ f4b7d3df-213c-46a6-8156-a83fcf6ffb34
begin
	
	# setup paths -
	_PATH_TO_ROOT = pwd()
	_PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
	_PATH_TO_DATA = joinpath(_PATH_TO_ROOT, "data")
	_PATH_TO_CONFIG = joinpath(_PATH_TO_ROOT, "config")
	
	# load required packages -
	import Pkg
	
	# If you do NOT have a Project.toml file -
	# Pkg.activate(mktempdir())
	# Pkg.add([
	# 	Pkg.PackageSpec(name="DataFrames"),
	# 	Pkg.PackageSpec(name="DelimitedFiles"),
	# 	Pkg.PackageSpec(name="CSV"),
	# 	Pkg.PackageSpec(name="Plots"),
	# 	Pkg.PackageSpec(name="TOML"),
	# ])
	
	# If you have already activated the project -
	Pkg.activate(_PATH_TO_ROOT);
	Pkg.instantiate()
	Pkg.update()
	
	using DataFrames
	using DelimitedFiles
	using CSV
	using Plots
	using TOML
	using Optim
end

# ╔═╡ 4b710f82-b1cc-11eb-2d14-11b6421a813f
html"""<style>
main {
    max-width: 1200px;
  	width: 70%;
	margin: auto;
}
"""

# ╔═╡ 0df967d7-ca1a-4a19-9c40-7f3a94764504
md""" #### Setup the Julia environment
"""

# ╔═╡ e0a279c1-7677-4c5b-a6e1-ee450b27020e
md""" ### Solution: Enzyme Activity Model

Derive and test an expression for the allosteric regulation of enzyme activity. 
Allosteric regulation is a fast mechanism that regulates the catalytic activity of metabolic enzymes.
In this type of regulation, the ability of an enzyme to catalyze a reaction depends upon the concentration of
effector molecules which do not directly participate in the chemical reaction. 

We will be developing a model for Phosphofructokinase (PFK) with and without the effector 3$^{\prime}$-5$^{\prime}$-AMP.

"""

# ╔═╡ aa6bba31-0d33-4476-acb8-707128e63a72
md"""##### a) Describe the Model
	
Let the model take the form $\hat{r}_{j} = r_{j}v\left(...\right)_{j}$
where $\hat{r}_{j}$, the overall rate of the PFK reaction ($\mu$M h$^{-1}$), 
is the product of a kinetic limit $r_{j}$ ($\mu$M h$^{-1}$), i.e., the maximum rate of conversion, 
and a control variable $0\leq v\left(...\right)_{j}\leq 1$ (dimensionless) that describes the influence
of effector molecules. Since we have only a single enzyme let j = 1. 

The model for the kinetic limit was given in the problem:

$r_{1} = k_{cat}E_{1}\left(\frac{F6P}{K_{F6P}+F6P}\right)\left(\frac{ATP}{K_{ATP}+ATP}\right)$

along with all the parameters associated with this model, and the concentration of enzyme ($E_{1}$), F6P and ATP. Thus, we can calculate the kinetic limit. However, we need to postulate a form for the v-varible, and estimate the parameters that appear in the v-variable expression. We'll do this from the data.

###### v-model
In the problem we proposed a three state model for $v\left(...\right)_{j}$: State 0: no effector+no substrate, State 1: no effector+substrate and State 2: effector+substrate. Following the procedure discussed in class, given this three state model, the v-variable takes the form:

$v_{1} = \frac{W_{1}+W_{2}f_{2}}{1+W_{1}+W_{2}f_{2}}$

where $W_{i}$ denotes the weight of state i (dimensionless), and $f_{2}$ denote the fraction of bound effector given by: $f_{2} = (x/K_{2})^{n_{2}}/(1+(x/K_{2})^{n_{2}})$. The quanity x denotes the concentration of the effector 3$^{\prime}$-5$^{\prime}$-AMP, $K_{2}$ denotes a binding constant (units of concentration) and $n_{2}$ denotes a binding cooperativity parameter (dimensionless). State 0 has no activity, while State 1 (with weight $W_{1}$) has basal enzyme activity, and State 2 describes the enhanced activity following from effector binding (with weight $W_{2}f_{2}$). 

 - The probability of state 0 is $p_{0}$ = $1/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.
 - The probability of state 1 is $p_{1}$ = $W_{1}/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.
 - The probability of state 2 is $p_{2}$ = $W_{2}f_{2}/Z$ where $Z$=$1+W_{1}+W_{2}f_{2}$.



"""

# ╔═╡ c4077a00-bef2-4003-ba51-8596550335b2
md""" ##### b) How do we get the parameters in the activity model?

We can estimate the parameters in the activity model from the experimental data given in the problem. This data set measured the rate of PFK at the same concentration(s) of ATP and F6P but different concentrations of 3$^{\prime}$-5$^{\prime}$-AMP. Thus, the kinetic limit was constant in this experiment. 
Thus, since the kinetic limit is constant, we can estimate a value of the v-variable as a function of 3$^{\prime}$-5$^{\prime}$-AMP:

$v_{1} \simeq \frac{\hat{r}_{1}^{*}}{r_{1}}$

where $\hat{r}_{1}^{*}$ denotes the _measured value_ of the overall rate. 

Lets look at the ends of the experimental data. Let $\alpha$ denote the ratio of rates when there is no effector, and let $\beta$ denote the ratio of rates when there is saturating levels of effector.

###### Case: no effector
When there is no effector, the v-variable model becomes:

$\alpha \simeq \frac{W_{1}}{1+W_{1}}$

which we can solve for $W_{1}$:

$W_{1} \simeq \frac{\alpha}{1-\alpha}$


###### Case: saturating effector

When the effector levels are saturating we expect $f_{2}\simeq{1}$ which gives the expression:

$\beta\simeq\frac{W_{1}+W_{2}}{1+W_{1}+W_{2}}$

which can be solved for $W_{2}$:

$W_{2} \simeq \left(1-\beta\right)^{-1}\Bigl\{\beta\left(1+W_{1}\right)-W_{1}\Bigr\}$

"""

# ╔═╡ 50886df5-d361-4646-898f-795a064c6001
md""" ###### Estimate the configuration weights
"""

# ╔═╡ a3264531-d9c8-42cf-b671-c6f46980a5a4
md""" ###### Estimate the binding and the cooperativity parameters

To estimate the remaining model parameters, we'll use an optimization approach e.g., Solver in xls or the [Optim](https://github.com/JuliaNLSolvers/Optim.jl) package in Julia. In this case, we want to minimize the squared difference between the model predicted v-variable and the v-variable estimated from the experimental data where the decision variables are $K_{2}$ and $n_{2}$:

$\min_{K_{2},n_{2}} \sum_{m=1}^{M}\left(v_{m,1} - v^{*}_{m,1}\right)^2$

where the summation is over the experimental data (the different values of the 3$^{\prime}$-5$^{\prime}$-AMP concentrations), $v_{m,1}$ denotes the model v-variable and $v^{*}_{m,1}$ denotes the experimentally estimated v-variable. 

Note: these estimates will __may not be unique__ because the problem is underspecified. Thus, your values may be different. 

"""

# ╔═╡ dc395064-55a7-42fe-8fb0-880692e7d81a
md""" ###### c) Plot model versus data. Does it fit?

Now that have the model parameters, let's plot the model v-variable versus the estimated v-variable and see if the model describes the data. In real life, we would do a statistical test on the residual (difference between the model and the data) and/or a cross-validation procedure (partition the data into test and validation sets) and calculate a training and prediction error. However, in this case let's just look - does the model generally describe the trend in the data?

"""

# ╔═╡ b1c395bd-6fa7-4e07-a339-894675791d91
md""" ##### Helper Functions Below Here
"""

# ╔═╡ 169f11a1-34dd-45a9-8f9a-6b18f4957998
function compute_activity_variable(effector_array::Array{Float64,1}, 	 parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    v_variable_array = Array{Float64,1}()

    try

        # process each effector -
        for effector_value in effector_array
            
            # compute -
            value = compute_activity_variable(effector_value, parameter_dictionary)

            # cache -
            push!(v_variable_array, value)
        end

        # return -
        return v_variable_array
    catch error
        throw(error)
    end
end

# ╔═╡ de9f574f-fc03-442c-90bd-810f44673466
function compute_activity_variable(effector::Float64, parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        W1 = parameter_dictionary["W1"]
        W2 = parameter_dictionary["W2"]
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # compute the f -
        f = ((effector / KD)^n) / (1 + (effector / KD)^n)

        # compute the v-variable -
        v_variable = (W1 + W2 * f) / (1 + W1 + W2 * f)

        # return -
        return v_variable
    catch error
        throw(error)
    end
end

# ╔═╡ 820c02ef-da91-482d-8bab-36de9b524fa1
function _obj_function(p,parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1}, 
    experimental_data_array::Array{Float64,2})

    parameter_dictionary["KD"] = p[1]
    parameter_dictionary["n"] = p[2]

    # compute the model v-variable -
    model_v_variable = compute_activity_variable(experimental_data_array[:,1], parameter_dictionary)

    # compute the residual between the data and model estimated v-variable -
    error = sum((v_variable_data .- model_v_variable).^2)

    # return -
    return error
end

# ╔═╡ b2fe595d-e90e-4041-bba4-8a1558bb83e2
function estimate_parameters(parameter_dictionary::Dict{String,Any}, v_variable_data::Array{Float64,1},
    experimental_data_array::Array{Float64,2})

    try

        # get unkown parameter guess from config -
        KD = parameter_dictionary["KD"]
        n = parameter_dictionary["n"]

        # setup the calculation -
        xinitial = [KD n]
        OF(p) = _obj_function(p,parameter_dictionary, v_variable_data, experimental_data_array)
    
        # call the optimizer -
        opt_result = optimize(OF,xinitial,BFGS())
    
        # return -
        return Optim.minimizer(opt_result)
    catch error
        throw(error)
    end
end

# ╔═╡ aa2267d4-ae47-4532-9fa0-460498b40fea
function load_experimental_data_set(path_to_data_file::String)::Array{Float64,2}

    try

        # load the experimental data -
        df = readdlm(path_to_data_file)

        # return -
        return df
    catch error
        throw(error)
    end
end

# ╔═╡ b25aeabb-2be1-48ef-bd72-3680ae8e9407
function load_configuration_dictionary(path_to_configuration_file)::Dict{String,Any}

    try

        # load the config TOML file -
        configuration_dictionary = TOML.parsefile(path_to_configuration_file)

        # return -
        return configuration_dictionary
    catch error
        throw(error)
    end
end

# ╔═╡ 5f846f30-f20f-4153-b87f-71ab429019cb
function compute_kinetic_limit(parameter_dictionary::Dict{String,Any})::Float64

    try

        # grab the parameter values from the dictionary -
        E1 = parameter_dictionary["PFK_concentration"]
        ATP = parameter_dictionary["ATP_concentration"]
        F6P = parameter_dictionary["F6P_concentration"]
        kcat = parameter_dictionary["kcat"]
        K_ATP = parameter_dictionary["K_ATP"]
        K_F6P = parameter_dictionary["K_F6P"]

        # compute the kinetic limit -
        kinetic_limit = kcat * E1 * (ATP / (K_ATP + ATP)) * (F6P / (K_F6P + F6P))

        # return -
        return kinetic_limit
    catch error
        throw(error)
    end
end

# ╔═╡ 5bab3b76-be8a-44c9-8de3-ac35992bda6a
begin
	
	# load the experimental_data file -
	path_to_experimental_data_file = joinpath(_PATH_TO_DATA, "Data-3-5-AMP.txt")
	experimental_data_array = load_experimental_data_set(path_to_experimental_data_file)
	
	# load my parameters -> configuration_dictionary contains, two sub-dictionaries 
	path_to_configuration_file = joinpath(_PATH_TO_CONFIG, "Configuration.toml")
	configuration_dictionary = load_configuration_dictionary(path_to_configuration_file)
	v_variable_parameter_dictionary = configuration_dictionary["vvariableparameters"]
	kinetic_limit_parameter_dictionary = configuration_dictionary["kineticparameters"]
	
	# lets compute the kinetic limit -
	kinetic_limit_value = compute_kinetic_limit(kinetic_limit_parameter_dictionary)

	# ok, so let's compute a data driven v-variable -
	data_driven_v_variable = Array{Float64,1}()
	(number_of_effector_values, number_of_cols) = size(experimental_data_array)
	for row_index = 1:number_of_effector_values
    	value = experimental_data_array[row_index,2]
    	scaled_value = value / kinetic_limit_value
    	push!(data_driven_v_variable, scaled_value)
	end
	
	# return -
	nothing
end

# ╔═╡ 2bb3a988-e054-4bfb-a306-cba7adb278f6
begin
	
	# compute W1 and W2 estimates -
	# W1 -
	⍺ = data_driven_v_variable[1]
	W1 = (⍺)/(1-⍺)

	# W2 -
	β = data_driven_v_variable[end]
	W2 = ((1-β)^-1)*(β*(1+W1)-W1)
	
	# return -
	nothing
end

# ╔═╡ 11df9c9f-e58f-450a-b630-3010db01391b
(W1,W2)

# ╔═╡ 8d8f4af4-8519-4a83-9027-bbf63b56e038
# estimate parameters -
parameter_array = estimate_parameters(v_variable_parameter_dictionary, data_driven_v_variable, experimental_data_array)

# ╔═╡ 81e97235-a5d6-4404-a701-c089f76aa791
begin
	
	# note: I put the values for K and n into the config file -
	
	# let's make an array of effector values -
	effector_concentration_array = range(0,stop=1.0, length=100) |> collect
	
	
	# compute the model v-variable -
	model_v_variable = compute_activity_variable(effector_concentration_array, v_variable_parameter_dictionary)

	# plot the data -
	plot(effector_concentration_array,model_v_variable,label="model",legend=:bottomright)
	scatter!(experimental_data_array[:,1],data_driven_v_variable,label="mean data")
	xlabel!("Effector 3-5-AMP (mM)", fontsize=16)
	ylabel!("Estimated v-variable (AU)", fontsize=16)
	
end

# ╔═╡ Cell order:
# ╟─4b710f82-b1cc-11eb-2d14-11b6421a813f
# ╟─0df967d7-ca1a-4a19-9c40-7f3a94764504
# ╠═f4b7d3df-213c-46a6-8156-a83fcf6ffb34
# ╟─e0a279c1-7677-4c5b-a6e1-ee450b27020e
# ╟─aa6bba31-0d33-4476-acb8-707128e63a72
# ╟─c4077a00-bef2-4003-ba51-8596550335b2
# ╟─50886df5-d361-4646-898f-795a064c6001
# ╠═5bab3b76-be8a-44c9-8de3-ac35992bda6a
# ╠═2bb3a988-e054-4bfb-a306-cba7adb278f6
# ╠═11df9c9f-e58f-450a-b630-3010db01391b
# ╟─a3264531-d9c8-42cf-b671-c6f46980a5a4
# ╠═8d8f4af4-8519-4a83-9027-bbf63b56e038
# ╟─dc395064-55a7-42fe-8fb0-880692e7d81a
# ╠═81e97235-a5d6-4404-a701-c089f76aa791
# ╟─b1c395bd-6fa7-4e07-a339-894675791d91
# ╠═820c02ef-da91-482d-8bab-36de9b524fa1
# ╠═169f11a1-34dd-45a9-8f9a-6b18f4957998
# ╠═de9f574f-fc03-442c-90bd-810f44673466
# ╠═b2fe595d-e90e-4041-bba4-8a1558bb83e2
# ╠═aa2267d4-ae47-4532-9fa0-460498b40fea
# ╠═b25aeabb-2be1-48ef-bd72-3680ae8e9407
# ╠═5f846f30-f20f-4153-b87f-71ab429019cb
