Modeling Psilocin Plasma Concentration Against Subjective Effects
Advisor: Movshovitz
Internship Location: Remote, CMU
Onsite Mentor: Hanyu Zhou, Graduate Student , Carnegie Mellon University
Hello! I'm Eileen Gong! My project aims to develop a mathematical model that predicts the intensity of psychedelic effects from psilocybin based off the concentration of its active component, psilocin, that is circulating in the bloodstream.
While current research has created charts detailing information of time against blood concentration and time against subjective effects, no one has created a comprehensive equation that connects the two individually. Using data from published clinical trials, I’m creating a pharmacokinetic-pharmacodynamic (PK-PD) model that traces the pathway, starting with the orally taken dose to blood concentration and finally to subjective effect
intensity. The main focus will be determining the time lag between psilocin levels peak in the blood and when effects peak in the brain. This research connects pharmacology with applied mathematics to create a framework for moving beyond dose-response relations to mechanistic, predicative modeling.
My final product will be a pharmacodynamic equation that relates oral dosage to blood concentration and a PK-PD model that relates blood concentration to subjective effects upon the brain. There will also be graphs showing model predictions compared to actual clinical data in order to demonstrate the model’s fit. Optimally, I will also create a dosage table showing “for average adult, dose x produces peak effects of y at time z” that researchers can reference for future psilocybin studies.
Project Blog
Post 1: Blog Post 1: Introduction
Welcome to my first blog post! In this one, I will be explaining the premise and goals of my senior project, along with my work in the first week.
Motivations
My Capstone class in differential equations was what led me to this project. I had done a class assignment on deriving the general Bateman equation and studying its uses. While researching, I stumbled upon how it’s applied in the field of pharmacokinetics, which I found incredibly interesting for how it treated the human body as a mere vessel through which substances traveled in and out of. So I decided to focus my senior project in this area, both to learn more about the applications of differential equations in pharmacology and to gain experience in creating my own models.
Psilocybin and Established Practices
Psilocybin is the compound found in numerous "magic mushrooms" that is known to cause hallucinogenic effects. Ongoing research focuses on studying its uses in therapeutic settings.
After psilocybin enters the body, usually through oral methods, it’s broken down into psilocin in the bloodstream. The compound directly affects serotonin receptors, most significantly the 5-HT2A receptor, leading to hallucinogenic effects. Psilocin plasma concentration peaks around the 1.5-2nd hour after dosing, and the volume of distribution (how far the psilocin travels beyond the blood) can go from 277 to 1016 liters, indicating significant tissue distribution.
Some effects of psilocybin include: distorted sights and sounds, altered depth perception, intense emotions, increased heart rate, nausea, inability to perceive time, increased openness to ideas, etc. Beyond just short-term effects however, it’s also been observed for various long-term effects that show up weeks after dosing, due to changes in the 5-HT2A receptor binding. This led to improvements in anxiety, depression, and other various mental disorders.
Many experiments have been conducted with psilocybin in varying doses. For Pharmacokinetic-Pharmacodynamic modeling (PK-PD modeling) approaches, it's been concluded that one-compartmental models (treating the problem as the drug entering and leaving a “compartment”) with first-order absorption and elimination kinetics (drug is absorbed/eliminated at a constant rate dependent on dosage amount) works best. The sigmoid Emax model is favored and is generally used in pharmacodynamics. It describes the effect of a drug for a specified concentration at the effect site (where the drug affects the brain), while taking into consideration the baseline effect (E0, the effects after one is given a placebo) and maximum effect (Emax, the maximum effect the drug can produce).
My Project
The question I'm researching is: What is the quantitative relationship between psilocin concentration in plasma and the subjective effects intensity following oral psilocybin administration?
To summarize, I'm aiming to create an equation that relates psilocin (the active component of psilocybin) plasma concentration to subjective effect intensity. Current research includes numerous published clinical trials detailing information on how psilocin blood concentration changes against time and how the intensity of subjective effects felt by the user changes against time. I hope to use this data to create a pharmacokinetic-pharmacodynamic (PK-PD) model. My final product will be an equation that relates the oral dosage to blood concentration and a PK-PD model that relates blood concentration to subjective effects upon the brain. There will also be graphs showing model predictions compared to actual clinical data in order to demonstrate the model’s fit. Optimally, I will be creating a dosage table that researchers can reference for future studies.
So what's the point? Practically speaking, this model could help researchers design psilocybin trials with optimized dosages. Although not immediately clinic-ready, it provides a framework that could reduce the 20-30% of patients who experience overwhelming effects without the added benefits. The equation I'm creating is the pharmacodynamic portion of the overall PK-PD system.
Week 1 Work
Over the past week, I have mainly focused on a preliminary literature review and building my foundational skills in Python for scientific computing. I've read four papers, all mainly centered around the results of past experiments and graphs/inferences created from that data. There are cases of independent experiments conducted by the researcher and cases of data being taken from past experiments. There was one paper that posed a similar reserach question as my work, but their result had an incredibly low accuracy. It serves as a good foundation and justification for the validity of my research. My main purpose in reading these papers was mainly to study how they performed their data analysis and also to collect the data necessary for my own anaalysis. I will be continuing this work and aim to finish litearture review by week 2. In terms of Python, I've been familiarizing myself with the NumPy, Pandas, and Matplotlib libraries by focusing specifically on the tools I'd need to perform my work. I will be working in Google Colab. I have so far ensured my foundational skills are solid by digitizing a concentration-time curve from one of the papers I've read and plotted it. Next week I will start focusing on the SciPy library, which is used for fitting curves.
Week 2: Literature, Data, and Learning the Tools
Hello everyone! Welcome to my second blog for this senior project, I’m so glad to see you all again!
This week, my main focus was on performing a thorough reading of existing literature to see how much of it I could use and the work required of me in order to analyze it, along with beginning to learn the computational tools I’ll need.
My original list of papers turned out to have very limited useful data, but luckily, I chanced upon a paper that had collected and referenced numerous other articles, and from that I gathered 8 experiments that met the criteria: oral psilocybin administration, healthy participants with reported ages and sample sizes, no substance abuse history, have both effect against time and concentration against time data, and most importantly, must have the participant’s subjective effects measured using the Visual Analog Scale (VAS). I was originally prepared to convert between subjective effect scales, but seeing how commonly used the VAS was I decided it was a good standardized scale to settle on. The only practical problem is that some of the papers don’t use the exact same VAS, as in some report in 0-10 while others report in 0-100%, so I’ll need to normalize this data to a common scale first.
Combined, these 8 papers provide me with data from 189 participants who received psilocybin doses ranging from 3mg to 30mg. I will be using all of this data, as my model maps psilocin concentration directly to the effect, independent of dose, and having a wide variety of data strengthens the curve fit. The low-dose subjects secure the bottom of the sigmoid curve, where effects are just emerging, while the high-dose subjects define the upper area. I mention sigmoid curve since that is what I predict the final model to take the form of, as drug effects typically follow the pattern of a slow onset at lower concentrations, an incredibly steep rise in the middle range, before a plateau at the higher concentrations.
On the technical side of things, I’ve begun learning two tools this week. WebPlotDigitzer is a browser-based tool that I’m hoping to rely on for data extraction. Most of the data I currently have is in the form of plots, and I will need to get it into a table format. I’ve practiced this so far on the data from Madsen et al. (2019) which show data from 8 subjects. Extracting the data from the two graphs and matching their times gave me my first real paired dataset. I’ve also begun working on tutorials for the SciPy library, mainly focusing on solve_ivp, which I believe will be at the core of my project. It solves differential equations, which is what a PK model is. I followed these tutorials and tested my knowledge by running the code in Google Colab to form simulated concentration-time curves and matching them with the data in published papers.
Next week, I will begin data extraction across the 8 papers and hopefully create the master dataset that will form the basis of my project. See you all then!
Week 3: Extracting Pharmacokinetic Data
Hi everyone! Welcome back to my third blog post. This week, I mainly focused on extracting the pharmacokinetic data necessary for my project. Allow me to break down for you guys what exactly that means and also what I’ve discovered along the way.
The primary tool I used for this, as mentioned before, is WebPlotDigitizer, a free web-based program that helps me extract data from a graph. You upload the image, calibrate the axes by clicking on known coordinate points, and then select specific datapoints which the program translates into real data. For graphs that showed only mean ± SD data, I ran three separate extractions per curve, one on the mean, one on the upper error bar, and one of the lower error bar, before calculating the standard deviations at each time point.
The only limitation I’ve noticed during this process, however, is the resulting data may not be incredibly accurate as I am manually selecting where I believe the center of the data point is, causing an inherent possibility for human error as a misclick of a few pixels is enough to change the extracted value. This is unavoidable since I’m working with image figures rather than the raw numerical data. Overall though, I’m not terribly concerned about this effecting my model, since the data representing 189 datapoints should average out across the entire dataset, rather than actually changing my results in a clear direction; as such the large sample size is incredibly necessary for this project.
One methodological change I’ve had to make was how my dataset is structured and the information it’d contain. Most of the papers I’m taking data from have only group mean data rather than individualized results, while some others do provide time data on each individual participant. I’ve decided to keep these two kinds of data separate, with the mean data being used to create my training set for fitting the pharmacokinetic model, and the individualized data will be kept out as a validation set. Meaning, model will never see the individualized data during training and I’ll use it afterwards to test how accurate the resulting fitted curve predicts individual-level observations. One thing I must address though is that this means I do not actually have data for 189 individual participants, but am instead working with groups means. The more accurate way to frame this is that my model will be the result of data representing 189 participants across 9 independent studies with doses ranging from 3mg to 30mg.
The data for concentrations is now extracted and organized. Next week I will begin the exact same process but for the pharmacodynamic data, meaning I will be extracting the data on the “any drug effect” graphs with VAS scores to form the other half of my master dataset. See you all then!
Week 4: Extracting Pharmacodynamic Data
Hey guys! Last week, I had just finished collecting the pharmacokinetic data. This week, I turned to the other half of my data, subjective effect intensity vs time. This is the pharmacodynamic side of my dataset, describing "how high people felt," that I'll eventually pair with the plasma concentration values in order to build my concentration-effect equation.
The extraction process was largely the same as last week, using WebPlotDigitizer to collect values from the 8 papers. However, this side of the data also came with its own set of complexities, each needing to be dealt with carefully.
The first was standardizing the intensity scale. Across my 8 papers, the experimenters didn’t always measure subjective effects with the same scale. Earlier I had mentioned that I looked for papers using a Visual Analog Scale (VAS), where participants mark a point on a 100mm line ranging from "nothing" to "extreme," producing a continuous score between 0 and 100%. There was one paper though, Madsen et al., that used a Likert scale instead, where participants chose a whole number between 0 and 10. These are not identical, but since they do use the same underlying construct and the conversion is straightforward, I stuck with it and just multiplied all the Likert scores by 10, ex. making a 7/10 into 70%. All subjective effect data in my master dataset is now on a 0–100% scale.
The second decision was which VAS measurement to use. Several papers reported many different kinds of subjective effects simultaneously, for example there was "any drug effect," "good drug effect," "ego dissolution," "drug liking," etc. I decided to exclusively use the "any drug effect" measurement across all papers. This is a very clear and simple measurement that captures the overall strength of the psychedelic experience, regardless of whether it was pleasant or unpleasant. It’s also convenient since it allows me to gather data from papers that didn’t classify their subjective effect intensity measurements as thoroughly. Using a specific subscale like "good drug effect" could have possibly introduced some bias, since an intense experience would be connected with high plasma concentrations but could also be categorized under low "good drug effect" ratings, interfering with the relationship I'm trying to model.
The same mean ± SD challenge from last week also applied here, with one addition. Some figures used graphs that had error bars leaving a measured value to represent confidence intervals, as opposed to standard deviations. A confidence interval indicates how precisely the mean was estimated, rather than the spread of the raw data. Since the confidence intervals in these figures were consistently narrow, though, that meant the mean values were generally well-estimated and reliable. For these figures, I extracted only the central measured value, while noting down that the error bars were present but small.
Again, the way WebPlotDigitizer is designed to be used manually does mean that small clicking errors are unavoidable, however, the same argument applies: with data representing 189 participants across 8 studies, random noise should average out rather than bias the results.
The master paired dataset is now coming to completion. Next week, I begin matching the concentration and effect time curves together and loading everything into Python for the first time. See you all then!
Week 5: Organizing the Dataset
Hey everyone! Back again. For week 5, my main focus was on organizing my raw extracted data into something the computer can process, something that can be read and curves can be made from.
My work was overall on reconciling the data points to a specific time point. My concentration data and subjective effects data were digitized from differing graphs within each paper, and said graphs did not always share the same time resolution. For example, a concentration measurement might exist at t=0.123 hours while the nearest effect measurement was at t=0.132 hours. These are close enough in time to represent the same moment biologically, but a computer would treat them as different time points and won’t pair them automatically. As such, my solution was to round all the time points to a consistent resolution, one whole number point in time ending with a 5 or 0, based off what seemed to be a reasonable pattern in the data collection. I did this manually, which was time consuming, but with the alternative being interpolating between time points mathematically, I took this route since the present data already had enough noise and this was more straightforward than adding another layer of estimation.
Once time points were reconciled, I classified each row in the master dataset with two specifications. One was the data type, either a mean of a group or from an individual. If it came from a mean I also recorded the group number so that when I fit the curves I can weight the data accordingly. The individual data is, again, reserved for my validation set. The second classification was whether a row had paired data, meaning it had both the concentration and subjective effect data for one time point. Rows with only the concentration data is still useful for pharmacokinetic modeling so they were kept.
The dataset also has a standard deviation column with every mean value. These SD values won’t be used for fitting directly, instead they will serve as an indicator for measurement uncertainty. I will re-run the curve fittings to the mean +SD and mean -SD bounds in order to understand how accurate my EC50 and Emax values based off the mean curve may be.
This week was just quite a lot of manual labor, but setting up this data set correctly is vital to my project. Next week, I will be starting work on creating and testing the pharmacokinetic model. See you all then!
Week 6: The Pharmacokinetic Model
Hey guys! This was definitely the most technically demanding step of the project so far, involving much more trouble shooting than I’d anticipated. However, after much trial and error, I’ve managed to create a working pharmacokinetic model for my data, aka a mathematical description of how psilocin plasma concentration changes over time after a specific dosage of psilocybin.
First off, I decided to derive some of the PK parameters myself instead of simply taking from the papers of my datasets. In this way, the model would reflect my pooled dataset much more accurately as opposed to just weighting and averaging established parameters, especially for a project focused on synthesizing data across studies.
The three main parameters I needed were the elimination rate constant (ke), apparent volume of distribution (Vd/F), and absorption rate constant (ka). I derived the ke by performing a log-linear regression on the elimination phase of my concentration-time data, or just when the concentration started to decrease. The natural log of concentration plotted against time produces a linear graph whose slope is the negative of the ke. I then derived Vd/F using the AUC (area under curve) trapezoidal method. I did this separately for each of my datasets, then took a weighted average of the results based on the number of participants for each study. I took my final Vd/F and ke to fit the ka for each dataset, using scipy’s curve_fit function, before taking a weighted average of the ka results. My final values are: ke = 0.3 hr^-1, Vd/F = 529 L, and ka = 0.714 hr^-1. The resulting model is a function that gets scaled by dose, as the ingested amount does affect the graph proportionally. I have graphed it against the actual datasets and it is a close match.
With the PK model now established, my focus shifts to the pharmacodynamic work, fitting a Emax sigmoid model that can map psilocin plasma concentration to subjective effect intensity. See you all then!
Week 7: Building the Pharmacodynamic Model
We’ve finally reached the main portion of my project, creating the pharmacodynamic (PD) model! Just to refresh, a pharmacodynamic model shows the relationship between the intensity of subjective effects felt by the user and the concentration of the active drug in their blood.
The model I used is called the sigmoidal Emax model, also known as the Hill equation. It describes the PD relationship as starting with minimal effects at low concentrations, before rising very steeply in the middle range, then plateauing at the maximum regardless of how much higher the concentration goes. The equation has three parameters: Emax, the highest felt intensity of subjective effects; EC50, the concentration at which half of that intensity if felt; and gamma, or the Hill coefficient, which determines how steep the middle portion of the curve is. The sigmoid Emax is the standard model for PD research, and looking at my aggregated dataset, the overall shape appeared to fit well.
I fitted my data to the curve by once again using Python’s scipy curve_fit function. The results were quite positive: Emax = 98.79%, EC50 = 11.22ng/ml, and gamma = 2.20. The model had an R^2 of 0.941 (meaning the model could explain 94.1% of the variance in subjective effect intensity) and an RMSE (root mean square error) of 6.82%.
One thing of interest to note is the gamma value of 2.20, indicating that the overall concentration-effect relationship is moderately switch-like. Simply put, effects are muted at lower concentrations, increase steeply after crossing a threshold, and then plateau. This is pharmacologically consistent with a drug acting through receptor saturation.
Next week, I’ll be looking at creating my final product: a dosage table. See you all then!
Week 8: Building the Dosage Table
Finally, we've reached the end product of my project, the result of all my hard work, the dosage table! To clarify, a dosage table is meant to tell you what amount of a specific drug you should take in order to achieve a specific subjective effect intensity.
The table was made through a two-step backwards calculation through the PK and PD models I've already built. The first step was to use the PD model to figure out the max plasma concentration for a specific target effect intensity. The sigmoid Emax can very easily be rearranged to solve for the concentration given the effect. I just plugged in my fitted parameters and solved. The second step converts the required concentration into the dosage through the PK model. However, rather than relying purely on the PK equation, I used a more direct approach by taking the observed peak concentration from each of my datasets. dividing that by the dose, then computing a weighted average. This gave me a scaling factor that represented how many ng/mL of psilocin each milligram of psilocybin produces at peak. This scale factor was a lot more reliable than a purely model-derived one, as it limits the exponential increase of errors, since it uses measured data as opposed to compounding on three estimated parameters.
The duration of how long the effect would last was calculated separately using the pharmacokinetic model. I put the threshold for "being high" at 6 ng/mL plasma concentration because that is what has been stated as the point at which users will feel the effects. Once the required dose is known, I simulate the full concentration-time curve for that dose and then measure how long the concentration stays above 6 ng/mL. This is a necessary step as it tells the clinician how long their patient will be experiencing the effects for and allows them to create an optimized therapy session. The resulting table is the following:
Target Effect (%) | Required Cmax (ng/mL) | Predicted Dose (mg) | Est. Duration (hr)
20 | 6.0 | 8.1 | —
30 | 7.7 | 10.3 | 1.6
40 | 9.4 | 12.6 | 3.1
50 | 11.3 | 15.2 | 4.1
60 | 13.7 | 18.3 | 5.0
70 | 16.8 | 22.5 | 5.8
80 | 21.7 | 29.0 | 6.8
90 | 32.3 | 43.1 | 8.2
The table overall does align quite closely with values reported in the literature. One thing to notice, however, is that the predicted dose for the 90% target effect seems to be unusually high. This is because the model is extrapolating beyond available data, as the highest dose in my datasets is 30 mg, and the highest reported effects were in the 80s. As such, there is not much available data for the 90% range, and so the model has to estimate based on the equations I created. Additionally, the 20% row shows an estimated duration of 0 hours, but that's mainly because its required Cmax itself is only 6.0 ng/mL, meaning there is likely very little time that the plasma concentration is actually above 6.0 ng/mL, and so the results say the estimated duration for a 20% effect is 0 hr.
The table is a population-level prediction, not a clinical recommendation. Individual pharmacokinetic variability means any individual's response could differ meaningfully from these predictions. What the table offers is a starting point derived from 68 participants across 4 independent studies, which is more than has previously existed in a single synthesized form.
Next week, I will be doing validation! See you all then!
Week 9: Validating the Results
Having the model and dosage table, I’ve come to the next most important step of this project, validating my results. Building the model is one thing, but demonstrating that it actually works is another, and arguably the most important scientific contribution. This week was dedicated to validation, specifically testing whether my PD model and dosage table held up when challenged against unknown data.
Unfortunately, I could not find another entirely independent dataset to perform my validations against. Instead, I performed something called “leave-one-study-out” cross validation, or LOO. The idea is simple: instead of fitting the model on all of my datasets, I take out one dataset, refit the model on the rest of the datasets, then validate that refit model against the dataset that was taken out. If the model performs well on data that it was never trained on, that is evidence of generalizability.
Overall, the results for the PD model were positive. Across all 4 holdout iterations, the model consistently reproduced the PD model with R^2 values close to the full mode’s value of 0.941. The EC50 and gamma parameters both stayed within a narrow range through all holdouts.
For validating the dosage table, I again used holdout validation, but instead compared the resulting dosage tables against the full model’s dosage table to see how far it differed. The results were quite positive, especially for target effects between 20%-80%, which had very small standard deviations, typically around 1-2mg. The 90% row showed the most variability across LOO iterations, with a standard deviation of nearly 5mg. However, this is again likely due to extrapolation, as during training the model only has data up to 80% subjective effect intensity to work with.
All in all, the validation results gave me reasonable confidence that both the PD model and dosage table reflect a reliable pattern in the data.
Next week I will be focusing on refining my project and creating my presentation,
Week 10: The Finale
Hey everyone! Don’t really have any substantial updates this week as I spent it mostly focused on creating and rehearsing my final presentation. The research portion is done, so now it’s time to take 9 weeks of work and shove it all into a 15 minute presentation. I organized the data, created slides, streamlined my process, and displayed my results. Wrote down a script, minimized it to within the time limit, and spent some time going over everything with my advisors, including performing a practice run.
For creating the slides, I decided I had to include a lot of visualizations in order to get my point across more clearly. PK and PD models are difficult to explain alone, and validation results are much more easy to comprehend the significance of in imagery. I created graphs and tables for my data. This included an image of my PD model model, overlaid with the actual datapoints color coded to match their datasets. Concentration runs along the x-axis and subjective effect intensity along the y-axis. I also created another graph that put the LOO curves in a light blue behind the main full model curve in red, to show the max and min of the hold out results.
The presentation is ready and the work can now be communicated clearly. I will do my best to walk the audience through my project and leave them with a good understanding about how psilocybin gets processed through the body.
Post 11: Afterword and Reflections
Hey guys! Just finished the “high” light of my project, my presentation! Overall, this project was a great experience for me, especially in terms of figuring out how to deal with data and seeing the methods that pharmacologists have developed to analysis drug processes through the body. It’s so interesting how what had started out as a scrapped idea from a differential equations class has now grown into three months of genuine research.
Of course, as all things do, there were the downs along with the ups. Hours spent combing through review papers trying to find usable datasets, the bouts of disappointment and annoyance upon finding a new public paper only to see reused graphs from data I already had, the effort spent in matching concentration and subjective effect datapoints by time and watching my data shrink as I removed unmatched time points. All in all though, the successes outweighed the losses. The amount of joy I got when I finally saw my numbers turn from “overdose” to “great experience”, when my graphs loaded and my model matched with the datasets, when my validations ran and the holdouts lined up with my full model, was just incomparable.
The presentation portion also helped me truly revel in success of my project. Communicating my months of hard work to my audience, watching as they listened intently, answering their curiosities, overall it was just such a satisfying experience. It is one thing to achieve success, it is another to fully understand your work and manage to minimize it down to a 15 minute, fully comprehensible presentation.
If I had the chance to continue with this research, I‘d love to incorporate individual-level raw data rather than just digitized published mean figures. That would greatly increase the accuracy of my PK model and perhaps give me more accurate dosage numbers that could actually be applied to clinical use.
Once again, I am thanking my advisors Dr. Movshovitz and Hanyu Zhou for supporting me throughout this project, for always being there for me and helping me through my struggles. And thank YOU all for accompanying on this journey, see you around!