How i can you ode extend function to used the solution till the end of timespan and in between used the RLC, offset and phase time as parameterizng function.

How i can you ode extend function to used the solution till the end of timespan and in between used the RLC, offset and phase time as parameterizng function. Kindly guide me.

[William Rose]

@Ehtisham,

Can you rephrase your question? I do not understand what you are trying to do. Of course you can increase the timespan if you want to extend the solution. I do not understand "in between used the RLC, offset and phase time as parameterizng function".

I have not run your code, but I have looked at it.

Others will be able to help you more, if you explain what your code is doing, and what you are modeling, in detail. A diagram would probably help. You are modeling receptor binding and unbinding, but I am not sure what else is going on. What are the "phases"?

Summary of what I infer from the code (wirttten after doing the analysis below): This code is a model of a receptor population. Each receptor is either active or inactive. The model predicts the number of acitve and inacitve receptors as funcitons of time. The conversion between active and inactive seems to follow standard first order kinetics. What makes it interesting and non-trivial is that the reaction rates are not constant. The reaction rates vary with time, and not in a simple way. Each reaction rate is constantly "relaxing" towards a target value, and the target values change for each phase of the simulaiton. There are 3 or more phases to the simulation.

In file (Run) file.txt:

You have

RLC = [0.1, 0.5, 1, 2.5];  % RLC Values 
% ...
% PhaseTimes should have one more element than RLC
PhaseTimes = [0, 100, 200, 300]; % phase times (each row defines a phase)

Why does RLC have 4 values? I expect 3, for R, L, and C. And, based on your comment, I expect Phase Times to have 5 elements, since RLC has 4. Are R and L and C circuit elements in a passive filter? RLC is not used in the differential equations passed to ode45().

PhaseTimes has only one row. Your comment suggests it should have multiple rows.

Please specify the units for all quantitites, including the rates. This can help with identifying errors in the equations.

The Tau values are all negative. The name "tau" suggests a time constant, and time constants are typically positive. Are the negative values correct? If you use the signs correctly elsewhere in the code , negative tau values are fine.

In file Code file.txt:

Function SimfileNPhase has 6 outputs and 13 inputs. It is a good idea to include comments at the start of any functIon, defining every input and output variable, including (a) what it represents, (b) variable type (scalar, numeric vector, string array, etc.), (c) units, if numeric. This will help you and others understand your code.

Function SimfileNPhase calls ode45(). I see that the state vector y has two elements: y(1)=Non_Active_Receptor_concentration and y(2)=Active_Receptor_Concentration. (It would be good to note this in the comments.) In many models involving receptors, the total number of each type of receptor is fixed, and the occupancy changes due to ligand binding and unbinding, or due to receptor conformational change. In such a situation, Unbound=Total-Bound, and you'd only need one differential equation for that part of the model. Do you really need two differential equations?

Function ode_LR() returns vector dy/dt. Examination of the code indicates that dy(1)/dt=-dy(2)/dt. Therefore the total number of receptors IS a constant, as I suspected above. This means you only need one differential equation, for either NonActive or for Active. Define a constant, TotalReceptors, in the main program, then call ode45() to solve for Active. When you get the results back from ode45(), compute NonActive=TotalReceptors-Active.

I am surprised that you wrote function findpeak(), instead of using the built-in function findpeaks(). I see that your function findpeak() finds the biggest positive OR negative peak, which findpeaks() doesn't do. Maybe that is why tyou wrote your own. But it is easy to find the largest positive or negaitve peak with findpeaks() (easier than writing your own version), and this will give you access to all the nice options which findpeaks() has. Using notation similar to yours:

[pkpos,locpos]=findpeaks(y,t,Npeaks=1);
[pkneg,locneg]=findpeaks(-y,t,NPeaks=1);
if pkpos>=pkneg
    Rpeak = pkpos;
    Tpeak = locpos;
    peak_type = 'High';
else
    Rpeak = -pkneg;
    Tpeak = locneg;
    peak_type = 'Low';
end

Function calculate_phase() has 4 input parameters, and one output, result, a structure with 8 elements. See comment above about explaining what each input and output is. I notice that calculate_phase() assigns the value RLC (which is one of the 4 inputs) to result.L. RLC is not altered inside the function, and is not used in any way, other than being passed back. What is the point of passing RLC in to the function, not using it, and passing it back?

I am sorry if my comments sound critical. I am trying to understand so I can help.

Answer 1

[William Rose]

@Ehtisham,

It appears to me that you are solving your system in successive time segments because you have parameters (RLC and offset) which change abruptly at certain times.

See the discussion here of how to extend a solution using ode45.

Instead of solving the system in spearate segments, you could apply the methods for solving odes with time dependent terms, discussed here. Perhaps you have already tried that, and perhaps you did not obtain acceptable results. In my epxerience, parameters that change instantaeously are a problem for ode45(), because ode45() reduces the time step around the boundary, in a futile attempt to find a smooth solution. Therefore you may want to use interp1() wth method 'pchip' to change the parameter quickly and smoothly around the transition time.

Kf and Kb decline exponentially toward their "target values". You could make Kb and Kf part of the state vector and compute them using a diferential equation, instead of using the nested functions which you now use. I think this approach would be more transparent.

Can you explain in words or with a diagram the role of offset, which is a part of the calculation of Kb? Why are there offsets for Kb and not for Kf?

[William Rose]

@Ehtisham, I will try. I will post more after I look at the code more.

[William Rose]

@Ehtisham,

[Edit: fix spelling errors]

I ran your code. I added the function in Code file.txt to the main program file and saved the resulting code as receptorDynamics.m. I ran the script and obtained a figure like the one you posted on Jan 21. I understand that you do not like the abrupt change in receptor concentration which occurs at the phase boundaries. In order to understand this better, we should examine the concentrations at the phase boundaries. Because you specify

timespan = 0:400;

you get reults at t=0,1,2,...,400, and not at intermediate times. It is useful to see the behavior in more detail. Therefore change the line above to

timespan = [0,400];

and rerun. Examine the results at phase change, t=200:

We see that the change in receptor concentration at t=200 is fast, but not instantaneous. It shows exponential dynamics. The time constant of the exponential is the time it takes for Concentration to go 63% of the way from the intial to the final value. Visual examination of the plot above indicates that this time constant is on the order of 0.01 s. We can get a more exact estimate of the time constant by

time constant=N/D

where D=slope of initial decay and N=difference between final and initial levels. We can do this as follows, after the script completes:

>> y=Active_Receptor_concentration;

>> N=mean(y(t>200.08 & t<200.1))-mean(y(t>199 & t<200));

>> tfit=t(t>=200 & t<200.004);

>> yfit=y(t>=200 & t<200.004);

>> X=[ones(size(tfit)),tfit];

>> b=X\yfit;

>> D=b(2);

>> tauFast=N/D

tauFast = 0.0102

This confirms our visual estimate of the time constant of the fast decay.

The model is behaving exactly as expected, based on the differential equations you have specified. During each phase, the model has two time constants: a fast time constant, tauFast=1/(Kf+Kb), and slow time constant which corresponds to the rate of adjustment of Kf and Kb. The plot of Kf shows the slow dynamics. Kb is similar.

The plot above shows that 100 s is not enough time for Kf to reach equilibrium, so its time constant must be on the order of 100 s or more. Perhaps you are surprised that the time constant for Kf is >=100 s. The code you use for Kf is

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - ...

Kf_LMax_values(i-1)) * exp(TauKf_ON * (t(phase_idx) - T_start));

You are using TauKf_ON as if it is a rate. Thefore the associated time constant for this equation is 1/TauKf_ON=1/0.01=100. If you want Tau_Kf_ON to be a time constant, then put it in the denominator, inside the exp() function, like this:

Kf_LMax(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - ...

Kf_LMax_values(i-1)) * exp((t(phase_idx) - T_start)/TauKf_ON);

and make similar changes for TauKf_OFF and for the code for Kb.

If you make the changes to Kb and Kf which I have suggested, and keep the paramater values unchanged, then the time constants for Kb and Kf will be similar to the time constant for Receptor_concentration. Threfore the dynamics of the two processes will no longer be separable, as they are in the current version. Then all process will be fast, and you will be able to reduce the duration of each phase to about 1 second or so, because this will be enough time for the system to reach a new equilibrium during each phase.

[William Rose]

@Ehtisham,

When I change the phase durations, I get unexpected results. For example, when I use the original phase durations,

PhaseTimes = [0, 100, 200, 300]; % phase times
timespan = [0,400];              % timespan

I get the results below.

Then I increase the phase durations to 500 s per phase, by changing the lines above as follows.

PhaseTimes = [0,500,1000,1500];    % phase times
timespan = [0,2000];               % timespan

Then I get the result below.

Examine Kb at the transition from phase 1 to phase 2, at t=100 or t=500. The plot below shows Kb from 5 s before the phase change to 55 s after. The top panel is for phase duration=100; the bottom panel is for phase duration=500. The vertical scale is the same for both panels.

In both cases, Kb=80 before the transition. In the first case (top panel), Kb goes down immediately, then rises gradually. In the second case (bottom panel), Kb goes upward immediately, then does not change. Also, note that Kb(t) exceeds the scalar Kb_max=80, for much of the time, in both cases above.

Is it what you want? If not, then find the problem and fix it. If you want help fixing it, then please explain in words the processes that govern Kb and Kf.

[William Rose]

@Ehtisham, I will help you, but I need you to explain the process that governs Kb and Kf. Do you have a journal article or another document that explains the system you are trying to implement?

[William Rose]

@Ehtisham,

Your explanation is very helpful. Kf_LMax depends on the receptor ligand concentration (RLC). The forward rate, Kf, should approach Kf_LMax exponentially, each time the RLC changes.

Please provide a similar explanation for the backward rate, Kb. In the code you provided, the asymptotic value of Kb in each phase is related to KbMax, KbMin, and offset, in a non-obvious way. The asymptotic value of Kb does not depend on RLC in your code. KbMax and KbMin are the same for all phases, but offset is different in each phase.

The bar plot below shows how KF_Lmax is related to the RLC values in the 4 phases.

Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 1, 2.5];  % RLC Values
Kf_LMax=Kf_Max*RLC./(1+RLC);
bar(string(RLC'),Kf_LMax);
xlabel('Receptor Ligand Concentration'); ylabel('KfLMax')

The screenshot below, which is part of an image I shared in one of my earlier comments, shows that Kf is not behaving the way you want. Specifically, at each phase change, Kf makes a very rapid jump followed by a gradualy adjustment. The rapid jump should not be there.

I will show you how to fix this, but first, please explain Kb, as you have done for Kf.

[William Rose]

In your recent comment, you have TauKf multiplied by time, inside the exponential. For example, you wrote

Kf_L(phase_idx) = Kf_LMax_values(i) - (Kf_LMax_values(i) - Kf_LMax_values(i-1)) * exp(TauKf_ON * (t(phase_idx) - T_start));

As I said in an earlier comment, I suspect this is wrong. Tau is usually a time constant. If your tau is a time constant, then you shoudl have time divided by TauKf, inside the exponential. But if you make the change, check the value of TauKf. You currently have TauKf=-0.01, so Kf adjusts with a time constant of (1/TauKf)=100 seconds. If you put TauKf in the denominator inside the exponential, then Kf will adjust with a time constant of TauKf=0.01 s, which is 10^4 times faster than in your simulations up to now. What is the proper rate at which Kb adjusts after each change in RLC?

[William Rose]

@Ehtisham,

This script shows a much simpler way to simulate your system. The state vector, y, has three components:

y=[ReceptorActive,Kf;Kb]

There are no complicated functions for Kf or Kb or Kb_Lmax, etc. In fact, the only function in the script is ode_LR().

When each phase is 100 s, this is the result:

When each phase is 500 s:

The unwanted jumps in the forward rates and the concentrations are fixed. The unexpected behavior when phase duration is changed is fixed.

You said "Same for the Kb just in the Kb we have a offset dynamic ranges", but your code treats Kb very differently than Kf. This is an insufficient explanation, so I have left Kb constant.

I suspect you want something like this for Kb(t):

where is the asymptotic value for the current phase, and depends on RLC (maybe) and offset for the current phase, and where may be different when going up versus going down. If that true, then specify the equation for . Because it's not obvious. In your code, the asymptotic values for Kb () do not depend on RLC, but the asymptotic values for Kf do. Is that what you want? In your code, the offset values cause an instantaneous change in Kb itself at each phase change - not just a change in the asymptotic value that Kb is going toward. Is that what you want?