Transit Compartment Model vs LAG Model (转移室模型与时滞模型)

有些口服药物在到达吸收部位前,在消化道内转运需要一段时间,这在PK曲线上,就会看到血药浓度在某某min以后才会开始上升。传统的处理方法(也是现在 主要药动学软件都内置的)是增加一个时滞参数TLAG。加上TLAG以后一般会得到较好的拟合结果,不过,从理论上分析,药物不能在TLAG以后,突然开 始吸收入血,肯定有个变化的过程。见下图的模拟,虚线是LAG模型,实线是TRANSIT模型:   转移室模型(TRANSIT)就是考虑了这个变化过程,用一些等体积、等转运速率常数的房室来逼近这个过程,原理见下图: 但是我们不知道存在多少个室(n)才合适,这需要通过模型拟合去求。若要用模型拟合,就需写出代数方程或者微分方程。可是,我们不知道有多少个室啊,也就是说微分方程不确定啊,貌似陷入死循环。但这些是难不倒大牛们的,看看他们是如何解决的。 其实只要知道最后一个室(an)中物质量的方程,即可把模型用微分方程表示出来了。经过复杂的推导,得出an的方程如下: 可是,这个方程是不连续的,难以用于求微分。这里作者采用了Stirling近似: 当n大于2时,Stirling近似与n的阶乘的差小于1%。近似以后的方程就是连续可导了,如下: 当然为了降低参数拟合的运算难度,可进行一下参数变换。 我觉得这种延时吸收的思路,不仅可以用于吸收过程的解释,也可以用在体内组织间转运,比如把传统的那种周边室,用一系列连续的转移室代替,转移室的最后又接回中央室,这样药物在兜了一圈以后又回来了,也许可以用来解释PK的双峰现象呢。 参考文献见附件: J Pharmacokinet Pharmacodyn. 2007 Oct;34(5):711-26. Epub 2007 Jul 26.   Implementation of a transit compartment model for describing drug…

一个用Gaussian计算溶剂中ECD光谱的方法来确定化合物绝对构型的实例

本文的例子是我们在研究鬼箭羽时发现的一个新化合物,当时其绝对构型采用经验规则与已知绝对构型的进行比较来确定的,发表在Natural Product Research上。后来,我们采用量化计算的方式,同样的确定了该化合物的绝对构型,下面详细的笔记,与诸位分享。Gaussian模拟溶剂下的ECD光谱 该化合物的英文名我们定为catechin lactone A (CLA),其化学结构和甲醇中实测CD光谱如下: CLA的结构和CD光谱(甲醇) CLA的2,3位为手性中心,通过它们在氢谱中的耦合常数,可知2,3位的两个氢处于对位,因此CLA可能为(2R,3S)和(2S,3R)两种绝对构型。下面详细说明具体的步骤: 1. 首先在gaussview中将构建CLA的3D分子模型,并对其结构进行优化。 这个化合物有两个手性中心,分别为图中箭头所指的位置。事先已经得知OH和Ar的朝向相反,所以若连接Ar的手性碳为S构型,则连接OH的手性碳为R,称为情况1.若连接Ar的手性碳为R构型,则连接OH的手性碳为S构型,称为情况2. 在情况1下,构建的输入文件如下: %nprocshared=16 %mem=32GB %chk=/home/xs/yjh/cdtry2-1C.chk # opt freq b3lyp/6-31g(d) scrf=(solvent=methanol) 0 1  C                 -5.55628238   -1.31201768    0.41502838  C                 -5.49337448   -2.68989041   …

基于pybel的在线化合物属性计算程序

在每天的裂解途径解析的时候,最后都要给出对应的化合物的结构信息,比如分子量、分子式等,每个化合物都要放到ChemDraw中去计算,然后再拷贝回来,效率低下,还比较容易出错。一直想写个在线的计算程序,输入SMILES,然后输出想要的结构性质和结构式,今天终于得空写好。思路如下:使用pybel的分子转化和属性计算功能,结合NCI的CIR和PubMed Compound,将这些综合起来,基本可以满足我们的要求。CIR的输入为SMILES,而PubMed的输入为InChiKey,所以需要用pybel将SMILES转化为InChiKey。具体代码如下: # -*- coding: utf-8 -*- import sys import string import cgi import pybel import urllib import openbabel as ob print "" print """ <!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML…

非线性拟合的一个前提条件:观测变量等方差分布

当我想好这个题目的时候,我顺便放狗搜了一下,找到几个非常对我胃口的博客(数据科学与R语言,Yihui Xie)。在看其中文档的过程中,我曾有短暂的失落感,“我在统计方面的知识太逊了”,不过我马上意识到自己又陷入了完美主义的坑,一个人怎能那方面都那么出色呢?在统计领域,我只要知道有哪些流行技术,哪些对我现在的研究课题有帮助,或者将来有帮助即可,只要分类整理入自己的笔记即可,完全可以等到需要的时候再去学。我自己喜欢的专业才是最根本的,这个是核心,其他只是点缀而已。我的目标是,“比药学家懂统计,比统计学家懂药学”! 为什么会想写这个题目呢?因为今天做酶动力学的模拟实验的时候,确实发现了这种现象,只有观测变量(因变量)的误差的方差恒定时,用非线性拟合得到的参数才是无偏的。对于线性拟合,只有满足高斯--马尔科夫定理条件时(零均值、等方差、不相关),最小二乘估计才是无偏估计。应该这对非线性拟合,同样适用吧。如果观测变量的误差不满足高-马条件时(比如方差随观测变量的增大而增大,可以从残差图上看到),该如何进行无偏估计呢?其实前几天已经提到了,正是Phoenix中的可进行群体药动学分析的PML!其所采用的技术为非线性混合效应模型法,广泛用于群体药动学参数的估计,模型中包含固定效应(相当于传统的非线性模型部分)、随机效应(个体间的变异)和误差模型(这里可以针对不同类型的误差进行建模,关键所在)。下面给出一个模拟实例,来说明用NLME进行非等方差时的参数估计。 模拟数据用EnzySolver的模拟功能来实现,具体参数如下: EnzySolver模拟数据,误差方差随浓度的增加而增加 模拟后的数据和误差分布如下: 用EnzySolver自带的非线性拟合功能,求得的参数为,与理论值有较大偏差,尤其是Km偏差快到10%了。将模拟的数据带入phoenix的PML模型,先假定误差模型为CObs=C+CEps (CObs观测值,C真实值,CEps误差项),拟合结果如下: 如上所求参数确有很大偏差,并且残差分布趋势随着浓度的增加而增大,这也提示所选择的误差模型不合适。将误差模型改为CObs=C+C*CEps,拟合结果如下: 可见,在我们选择了合适的误差模型以后,所求算的Km和Vm与理论值就非常接近了,且所估计的方差0.013与理论方差0.02相当接近。 –EOF– 文章来自,微信号:MS4Fun,不定期发布自己在质谱应用和建模&模拟方面遇到的一些有趣的事情,欢迎分享与推荐。

酶动力学的PML模型

群体的加权残差与观测浓度的分布 前面我们用Excel和Matlab实现了酶动力学的参数估算方法,今天我们介绍使用药动建模工具phoenix来进行酶动力学分析。正如我们前面提到的,phoenix不支持隐式方程的求参,所以需要对数据进行一下变通,通过观察IMM方程,我们发现,虽然C0和C1之间不能用显示方程表示,但是C1和t之间的关系却可以写成显示方程,一系列的C1和对应的反应时间t我们都知道,C0可以类比于PK中的给药剂量。这种处理方式类似群体PK,每个个体在t时间采一个点,存在大量个体,每个个体的给药剂量不同,这种稀疏的数据只能用群体的模型去估计。为了实现群体PK,这就需要用到PML,即Phoenix Model Language。PML是借鉴了C++,S-Plus中的一些规范的用于Phoenix平台的建模语言,可进行非线性混合效应(NLME)分析。NLME背后的原理理解起来有点困难,但我们作为药学人士,没必要搞懂其具体的原理,只要我们会用这门语言编写满足我们要求的模型代码即可。我的原则是,够用即可,多了不学。PML支持微分方程和显示方程,C1与t的关系用微分方程表示更简单,其实就是最常见的米曼方程: dC/dt=-Vm*C/(Km+C) EnzymeKinetics(){ //米曼方程 deriv(C = - Vm*C/(Km + C)) //dose的位置,即初始浓度是加在C这个变量上 dosepoint(C) //误差模型,一般是additive,看看残差分布是否随机,如果不随机,可以用multiplictive等其他误差模型 error(CEps = 1) observe(CObs = C + CEps) //CObs是观测值,这个不能变,其等于C加误差 /*结构参数,这个主要是用来考虑个体之间的波动,其认为存在一个群体的值,即tvVm和tvKm,而对于每个个体来说,它们又有自己特定的值*/ stparm(Vm = (tvVm) *…

隐式方程的matlab实现:酶动力学(Integrated Michaelis-Menten Equation)的无偏参数求解

神奇而强大的计算分析平台 Google搜索了一下,发现隐式方程的求参问题有不少解决方案,matlab, R都可以实现。基本思路与我用Excel求参的方法一样,先编写一个函数来实现隐式方程的显示化,matlab中要用到fzero或fsolve(Excel中用的是单变量求解插件),然后再用非线性求解函数进行参数拟合。下面是积分后的米曼方程(IMM)的matlab求参代码: 先构造IMM函数: function y=IMM(p,x) %输入参数 Vm=p(1); Km=p(2); y=zeros(size(x)); NN=length(x); %不输出提示 opt = optimset('display','off'); for i=1:NN %就是用这个方程来求出某个x(i)下的解,y的初值设定为x(i)*0.9 y(i)=fzero(@(y) x(i) -y-Vm*15+Km*log(x(i)/y), x(i)*0.9, opt); end end 上面的函数写好后,保存,最后只要用下面一行命令就可以求解参数了: X=lsqcurvefit(@(params, xdata) IMM(params, xdata),...…

酶动力学数据分析工具:EnzySolver 1.0

EnzySolver 1.0界面 前段时间进行体内外相关性建模,遇到需用体外酶动力学参数预测体内的消除程度,在与做试验的同学讨论的时候,发现大家应用米曼方程并不注意前提条件,比如需控制反应时间或起始浓度,使反应结束时,底物的消除不要超过10%,否则计算的反应速率就与真正的起始速率有很大偏差,用这样获得的有偏差的起始速率与起始浓度求算Km和Vmax,必定难以获得真实的酶动力学参数。 从纯数学的角度讲,这种使用最多的用起始速率和起始浓度来求Km和Vman的方法,无论反应时间多少都存在偏差,只是反应时间短的时候偏差小而已,我们把这种方法本身引起的误差成为系统误差。那是不是存在一种没有系统误差的方法呢?关于酶动力学参数的计算方法在上世纪70-90年代,研究最多,但那时计算技术并不发达,所以研究者们设计了多种米曼方程的变体,还有很多作图的方法来求Km和Vmax。我认为现在的计算技术已经足够发达,肯定能实现没有系统偏差的计算方法。还有一种酶动力学研究方法,就是测定随着反应时间的延长,底物浓度的变化,再通过拟合积分以后的米曼方程进行求参。但随着反应时间的增长,酶的活性会发生变化,所以这样测定的浓度,越到后面越不能反应酶的活性。从实验的角度考虑,还是测定系列起始浓度,反应固定时间以后,测定剩余浓度(或代谢物浓度)的方法更靠谱。 假定起始浓度为C0,反应t时间以后的浓度变为C1,改变起始浓度C0,我们就可以得到系列的C1。根据积分以后的米曼方程,可以得到C0与C1的关系如下: Integrated Michaelis-Menten Equation (IMM) 这个方程为隐式方程,不能表示为显示方程,也就是y=f(x)这种形式。现有的建模工具,就我所知,都不支持隐式方程的求参,也许存在这样的工具,但我不知道。我直觉上感觉用Excel可以实现,所以在一个月前开始用Excel的Solver和单变量求解工具设计能够对这样的隐式方程求参的算法,大约用了1周多时间,就写好了一个算法。模拟的结果显示,当数据中不含噪音时,可以无偏的求到正确的参数,而一但增加噪音,求参结果就显得非常不可靠,我一度认为这正式IMM的缺点。此外,我还设计另一种与现有的算法(TM)偏差相反的算法PM,希望通过两种算法的平均,达到偏差消除,从而求出正确的解。就在昨天,我还设计了Robustness参数,APDD, APF等参数,通过比较它们之间的大小,自动判断哪种算法获得的结果更可靠。本打算今天下午的Seminar向老板汇报工作,但老板有事外出,Seminar被取消。早知Seminar被取消的话,我这两天就不用那么辛苦赶工了,心里有种极放松的感觉。暂时没有其他的事情,我就在EnzySolver上模拟不同的数据,看看结果如何,可能存在什么规律。突然,我发现前面设计的求解算法有问题,当C1这一列没有噪音的时候,只要一个循环就可以得到正确的解,而当它有噪音时,整个算法Solver仅相当于只迭代一次,怪不得误差会比较大呢。在晚上回家前,终于完成了那个迭代算法,模拟结果显示,即便存在很大噪音,其求的参数总是在各种算法中最接近真实值的,竟然比我用群体的方法(理论上考虑了误差模型)还精确,不可思议。我本来打算用Phoenix的NLME方法作为压轴的终极解决方案呢。我暂不知道我设计这种隐式方程的求参方法是否别人已经发表过,也不知道它是否可以用于其他隐式方程,是否具有通用性,这些都需要后面的工具进一步验证。 –EOF– 文章来自,微信号:MS4Fun,不定期发布自己在质谱应用和建模&模拟方面遇到的一些有趣的事情,欢迎分享与推荐。

你就是自己的建模师

Snoopy是谁的模型?模型是对事物本质的逼近,有时我们宁愿选择不那么逼近真实,因为这样我们才能更不受世俗约束的表达自己的真实本质。 ​        在高楼大厦林立的繁华城市中,一只小狗躺在自己的屋顶上,仰望天空,陪伴它的是一只忘记了季节的小鸟“糊涂塔克”。它幻想开飞机,开火箭,当超人,它从不担心是否需要考飞机驾照,是否要学航天知识。。。。。。这就是Snoopy! ​        我不知道为什么会选择药学这个专业,也不知道为什么读了硕士、博士、博士后。生活所迫?也许吧。那时候读书是唯一能让我生活不太依赖家人的方法,不用交 学费,还有生活补助。但内心那种冲动这么多年来一直与现实抗争着:本科时,我着迷QBasic的图形编程,泡在东大的计算机中心,忘记了时间;硕士时,我 迷上了汇编语言,抱着看雪的加密与解密在不能站直身体的春运列车上看得忘记一切;博士时,我学习了matlab, C++, perl,并用perl开发了一个专业网站,稍有盈利;博后时,我迷上了建模,把WinNonLin用的纯熟(有点吹牛啊),但我很少用它解决药学内的问 题(不务正业),如果有人让我帮忙分析药动数据,我建议他们还是去用3P87或DAS之类的软件去算吧。但是如果你问的是WinNonLin从没有人用过 的应用,我会比较感兴趣, 比如有个博士生让我模拟几个中药成分的化学反应动力学,或者如何解决有损失采样造成的系统偏差,或者养殖鱼苗中给药后药物多长时间可以清除等等,我会比较感兴趣。这就是有点叛逆的我吧! ​        昨天我把放下好久的NONMEM,WFN,XPOSE又配置了起来,因为这就是我现在的工作!这是一个结束,也是一个新的开始。 ​        分享一下自己在08年研究NONMEN的时候写的笔记,正因为有了这个,我很快就恢复了建模环境: ================我是可耐的分割线================ 用G95编译NONMEM6 ​        NONMEM6的源码可以用Compaq Visual Fortran 6(CVF)直接编译,没有问题,若是用其他编译器,比如G95,就需要对源代码稍微改动一下,因为不同的编译器对一些命令的使用有区别。 ​        G95是免费的Fortran编译器,安装后只有十几兆,在网上很易找到并下载。而CVF的安装文件近600兆,很难在网上找到下载,并且安装后,占硬盘空间也很大,至少200多兆。因此使用精巧、免费的G95编译器进行编译具有很实际的意义。 ​        整个编译的过程可以参考NONMEM6的CDSETUP6.BAT文件,至少稍微修改就可以对NONMEM目录中的*.for进行编译,G95编译后的文件后缀为*.o,而CVF编译的后缀为*.obj。我修改后的setupG95.bat如下(注意,先删除所有目录中的*.obj,以及nonmem.lib,并删除run目录): 1. 运行setupG95.bat前,需要将根目录下MAKEG95BAT中的第45和48行中的“.obj”改为“.o”; 2. NM\FLU.for中的第9行“CALL COMMITQQ(I)”,改为“CALL    FLUSH(I)”,否则编译后运行会出现错误“c:\nmvi1\NM\nonmem.lib(FLU.o):FLU.for:(.text+0xa): undefined reference to `commitqq_'    ”; 3. NM\BLKDAT.for中的第28行中的'con'改为'CONOUT$ ';…

Phoenix实现“男生追女生的数学模型”

由于当前男女比例失调,男生能否追到女生具有重要现实意义。此模型有助于男生理性把握时机。 Phoenix是Pharsigh开发的药动建模平台,是制药界的“Gold Standard”。Phoenix不仅可以研究药动学问题,对于“男追女”的问题研究也是神器一枚。且看我们如何用它来实现获得了2013年菠萝科学奖的论文(“男生追女生的数学模型”)中的微分模型。 首先建立一个新Project “Boys after girls”,然后在Workflow中新建一个"WNL5 Classic Modeling-->WNL5 ASCII Format"。这是与旧版本兼容的建模语言,Phoenix还有自己更强大的PML建模语言,这个我们后面再讨论。此外,在Data项下,建立一个包含“Function, Time, Value”3列的表,function有2个,分别对应论文中的学业成绩函数y(t)和疏远度函数x(t),Time是1到500等间隔的数值,步长为1,步长越小,模拟的结果越精细,Value空着即可,是我们要通过微分方程模拟的值。 在窗口右侧的Setup中,写入下面的ASCII模型,其实绝大部分是套话,只有DZ(1)和DZ(2)两行代码最重要,对应论文中的方程组(5),我用Vw代替了“V1+lw”。 在Initial Estimate中设定f1,a1,k1,k2,m,n,Vw的初值(f1=0.03, a1=0.07, k1=0.954, k2=1, m=0.4, n=0.48, Vw=0.48, 在论文中有),最后勾选Engine Setting中的“Simulate”。 至此,男追女的数学模型建立完毕,运行后就可生成论文中图1-图4的数据。下面是图1的结果: 从图中可以看出,只要自己的能力和时间把握的比较好,就可以实现学习和爱情的双丰收(两个相互促进,共同提高)。祝看到此文的单身男士都找到自己心仪的女生! --EOF-- 文章来自,微信号:MS4Fun,不定期发布自己在质谱应用和建模&模拟方面遇到的一些有趣的事情,欢迎分享与推荐。