I am a relatively new user to Mathematica. Could someone tell me the best and most efficient way to code in Mathematica these three differential equations, shown in the image, for a three gene negative feedback loop? I have included an image below of how I was doing this previously with brute force in Excel using small time increments of 0.000025 min. Obviously, Mathematica is a much better way to do this calculation. You can see that each differential equation includes the gene product calculated by one of the other two differential equations so all 3 equations are interlinked. I would like to be able to simulate various three-gene negative feedback systems by changing the parameters on the left-hand side of the image below.
Thank you in advance for any help that you can provide me on coding in these three differential equations into Mathematica and informing me on the best way to do simulations in Mathematica.
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.
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.
I have done some more research on how to simulate with three coupled differential equations, but the code i have written in the attached Mathematica notebook is still not working. I am open to any advice to successfully simulate these three differential equations in Mathematica.
Then I see your three plots. But the plots do not look very close to the plots you show from Excel. I am guessing that this may be because of the values for some of your coefficients when you tried to translate from discrete steps in Excel to derivatives in Mathematica. Perhaps you could estimate the derivatives by looking at the change in values from 0.0 to 0.000025 and the slope of those changes might be close enough to the derivatives that you need to give NDSolve
Thank you immensely for the help on this! Much appreciated. I did have a few of the parameters not set correctly like K2=5e7 not 1e7 and the starting condition for x[0] should have been 1e-6 not 0. However, the bad news is that even after setting those parameters and initial conditions correctly and setting the range from 0 to 0.000025 the curves do not match the Excel curves at all as you pointed out.
I think your suggestion is a good one if I understand it correctly. Do you mean NDSolve in Mathematica can be given a time increment like 0.000025 to do the numerical solution with or are you suggesting something else? At this point, I am not sure I am capable of doing either of your two potential solutions because I am not an expert in Mathematica.
Could you provide a little more detail on any ideas you have of overcoming the differences seen in the curves between Mathematica and Excel? It is hard for me to believe Excel would be better than Mathematica at this calculation, and it is more likely my incompetence in Mathematica may be causing the differences we are seeing in the curves. I am pretty confident that Excel (while taking forever with one million rows to do the calculation) is producing the correct curves.
What is most puzzling is that your graph appears to show that two of the protein concentrations approach the same steady state value. It is difficult to see a justification for that in the equations. Is it possible your spreadsheet does not correctly implement the differential equations shown?
Beyond that, I suggest that if Mathematica results (using adaptive time steps) differ from those in Excel using a fixed time step, then Mathematica is more trust worthy.
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.
While x and z do have the same lifetimes, the coupling of equations breaks the symmetry. x’[t] depends on z[t], but z’[t] depends an y[t]. The equations with the parameters substituted look like this:
If you set the derivatives equal to zero, and solve for the three concentrations, they are not the same. (That is what DFixedPoint does. Its new in 14.1.) Maybe I have made an error in the equations or parameters.
Regarding time stepping, I did not specify time step parameters to Mathematica. It has very clever algorithms for doing that.
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.
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.
The equations in my notebook are the differential equations that I took from a textbook, Molecular Biology of the Cell, a classic in the field. I choose all parameter values myself (because the textbook does not give any particular values in the figures they show) by trial and error to generate an oscillatory behavior. It is well known that this negative feedback gene system described by the three differential equations should give oscillations in protein concentrations until steady-state is reached.
This seems like too easy of an answer to your question, so I may have missed what you are really asking me. If so, please get back to me.
Due to my inexperience with Mathematica, the last notebook I sent you did not evaluate correctly. In the newly attached notebook, i do see all proteins overlapping with one another, so that is pretty convincing evidence to me that Mathematica is working correctly.
i wish i was a better mathematician because I still do not understand why Excel is not giving the correct curves. The only think I can think of is that 0.000025 is still not a small enough time increment to allow Excel to simulate the differential equations accurately.
In any event, I think I will believe the Mathematica curves more than the curves in Excel.
Thank you again!
I think you had some previously defined definitions still in the kernel. That can be deceptive. If you quit the kernel and then evaluate the notebook it will be done with a fresh kernel and you should see the same results.
Any chance you could scan the page with the equations from the text and attach to a post?
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.
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.
I did misunderstand what you meant then. You were asking about the equations I used in Excel. Sorry that I totally missed that. I have attached the section below from the textbook dealing with these oscillatory, negative feedback systems. I did make some comments in my response to another user, David Keith, about that section, see Figure 8-80.
It will be interesting to see what the small time increment, step-wise Excel equations return in Mathematica. In Excel it takes about one million rows to get 30 min. of data with 0.000025 time increments.
Thank you for the curiosity and taking a look at those equations
