David,

Thank you for all the help. I will take a look at the notebook you attached. Your analysis is much appreciated. To answer your question about two of the proteins reaching the same steady-state--that is correct because two of the proteins (x and z) have all the same parameters so upon reaching equilibrium (steady-state in essence) their protein concentrations should be the same. For some reason, protein x and z are not reaching the same steady-state level in Mathematica. In regards to your second question, I did find an error in my Excel sheet (I attached a truncated version yesterday, which has the mistake in it). However, the mistake did not change the oscillatory nature of the curves and the correct curves are almost identical those you saw in my post 2 days ago, so I am still perplexed why I am not seeing the same behavior in Mathematica. The curves should oscillate and two of the proteins should reach the same steady state. I have attached a corrected truncated version of my Excel sheet in case you care to look at it. I can find no mistake in the way I did the analysis in Excel other than it is not eloquent and extremely slow due to 1 million rolls in Excel. If I understand what you did in your notebook, I believe you just told Mathematica how large of a time increment to take when solving the coupled differential equations; is that correct? If that is correct, it is even more puzzling to me because the increments appear very small, which is similar to what I used for time (0.000025 time increments) in Excel.

Scott

Attachments:

Oscillating equa...xlsx

David,

Thank you for the input. It makes sense what you are saying. I went ahead and made all the parameters the same for each type of parameter for each of the three proteins. Under those conditions, all the proteins should reach the same steady-state level. You will see in the attached notebook they do not, so I am still puzzled by what is going on in Mathematica. When I simulated in Excel, all three proteins reached the same-steady level, which should be the case since all proteins are being made and broken down at the same rate. I wanted to do this comparison to try to figure out why there is disagreement between Excel and Mathematica.

Any further feedback would be appreciated.

Scott

Attachments:

677785769-DK SK ...nb

David,

yes, I screwed up the evaluation of that notebook--our messages crossed on the wire. It will not be the first and last time I do that I am sure! Mathematica is a real insider's software in my opinion.

I have attached the pages you requested, see Fig. 8-80. You can forget the hill coefficient exponents in the text. Hill exponents are for systems with positive or negative cooperativity (where binding of one ligand either increases the affinity or decreases the affinity of the next ligand binding event). This current system is complicated enough without throwing Hill exponents into the mix. That is easy to do later once the foundation is laid, which I think we are finally there now from the help of users on in this community like you. I am very grateful for the help on this.

Scott

Attachments:

Molecular Biolog...pdf

Oops -- you were typing this while I was responding above. ;-)

Hi Scott, I’ll have a look at this tomorrow.

David

Thank you for all your expertise and help! I am glad it is fun for you. Not so much for me yet. Hopefully, someday I will get there. I will take a look at the notebook in an hour or two. I am really glad you engaged in this help request.

Scott

Do you have the original differential equations and parameters that you used to develop your discrete model? That might make it easier to produce the differential equations that MMA can accept, solve and plot. Or did you develop the discrete model from scratch?

Both Excel and MMA can be used to approximate things like this. The work is trying to translate the information into a form that each of them can use.

I've tried and can't get my hands on a copy of your textbook at the moment.

If you could copy the differential equations from the text into a message here and include the values of the parameters that you are using at the moment then I will try to translate that system into Mathematica and plot the solutions and see how close my plots are to your posted Excel plots. I might be able to generate your entire Excel sheet from the example you gave, but I'm not sure I can do that correctly.

Hopefully we are going to get you what you need. Thanks for your patience and assistance with this.

Bill,

Thank you for all the work! I chose Kx=1e8, Ky=5e7, and Kz=1e8, which could be realistic affinities (1/K) of 10 nanomolar and 20 nanomolar. Figure 8-83 is for the two-gene protein negative feedback loop. i chose the three-gene negative feedback system because the more genes the more the system oscillates. A five gene loop would have been even better, but I did not want to get carried away. You can leave the Hill coefficients (hx or y) out, if you would like, because a system with no cooperativity has hx or y=1. This system is complicated enough without more complexity. I choose all those 12 parameters arbitrarily though trial and error to create an oscillation the best I could in Excel. It will be very informative to see what happens in Mathematica when you enter your equations.

Scott

With that fine a time step and reasonably well behaved equations I doubt that the problem is in the accuracy of Euler's method.

Hi Scott,

I implemented Euler's method in Mathematica using a function called NestList. It repeatedly applies a function to the last result, starting with some initial condition. I used your equations for the derivatives and your time step. The original result that both Bill and I have using NDSolve is this:

The result using Euler's method is this:

You can see that they are quite similar. And they converge on the same concentrations. I suspect the difference arises from the fixed time steps being to large in the early part of the simulation, resulting in under damping. Mathematica has adaptive time stepping to prevent that.

I don't know why the Excel Euler's method produces different results, but I would trust NDSolve over the perhaps-too-simple Euler's method.

Attachments:

By Euler.nb

David,

This is beautiful work and what I have been after for about 5 days now! To me it seems clear from your two simulations using Euler's method that one never knows when the time increment is small enough to ensure accuracy. One also runs into computational time problems in Excel with really, really small time increments. I guess I would have to go to 10 million rows in Excel to get the accuracy I need. Thank you again for all the time you put into solving this conundrum for me.

Do you mind sending me your notebook where you implemented Euler's method? I would like to know how to do that method in Mathematica for comparison purposes.

Scott

David,

Sorry that I missed the attachment.

This whole issue is exactly why I tried to do this in Mathematica to see if there were differences and there are. I wanted the most accurate answer I could get and Mathematica gives me that.

That Euler notebook in Mathematica looks complicated. I have some studying to do there!

I am satisfied my original question has been answered thanks to both Bill and you!

Scott

David and Bill,

Here is a Youtube link showing the method (Euler's) I used to solve those three differential equations in Excel. https://www.youtube.com/watch?v=vXWb1Kspx8s I think you surmised this already, Bill, when you deconstructed my Excel equations. Euler's method is crude and slow, but I always thought if the dt was small enough, it approximated derivatives pretty accurately (hence my small time increment of 0.000025 min. and my 1 million row Excel sheet).

Scott

Bill,

I have attached a truncated version of my Excel file (the complete original file is 129 MB with 1e6 rows in it). You can see the formulas in the sheet that I used to simulate the concentrations of proteins x, y, and z from the differential equations for x, y, and z formation with respect to time. x, y, and z are the proteins made in the negative feedback system from the three genes, x, y, and z. All the parameters are defined in the Excel sheet too that are used in the three differential equations. Thank you for any help you can render to me.

Scott

Attachments:

Oscillating equa...xlsx

x,y, and z are the three different proteins that are made by the three-gene negative feedback loop and described by the three differential equations. Thanks for any help you can give me.