Dette projekt lavede jeg da jeg gik i gymnasiet. Projektet er kraftigt inspireret af en øveeksamensopgave vi lavede på Biologi A om hormoner.
Opgaven handlede om de tre hormoner LH, FSH og TSH, som der er en interessant sammenhæng mellem og som dette projekt handler om.
Af hensyn til den nysgerrige læser, vedhæfter jeg hele opgaven herunder. Det er dog kun spørgsmål 4 der er interessant i denne kontekst. Det har jeg sat en rød ring om, så resten er strengt taget ikke nødvendigt at læse for at forstå mit projekt, men kan gøre det lettere at følge med:
Opgaven går ud på at forklare resultaterne af et eksperiment. I eksperimentet indsprøjtes 100 mikrogram af et hormon (LH-releasing-hormon), hvorefter koncentrationen af et andet hormon (LH) måles løbende. Målingerne er plottet i et koordinatsystem vist i Figur 3.
Opgaven er interessant af flere årsager. Det er selvfølgelig spændende hvordan kroppen fungerer, men det der fascinerer mig er muligheden for at løse opgaven med matematik og programmering. For at se, hvorfor det er smart og interessant at bruge matematik til opgaven, vil jeg forsøge at løse opgaven uden brug af matematik.
Der er flere måder at forklare resultaterne af eksperimentet på, men det der sker er grundlæggende et eksempel på negativ feedback.
Det der skal forklares er, hvorfor LH-koncentrationen ændrer sig. Der ses en kraftig stigning i koncentrationen af LH kort efter indsprøjtningen af det andet hormon (LH-releasing-hormon), og senere ses et fald i koncentationen af LH. Forklaringen på stigningen i starten er, at det indsprøjtede hormon (LH-releasing-hormon) har en stimulerende effekt på dannelsen af LH, se opgavens Figur 2. LH-koncentrationen går ned igen efter lidt tid. Det skyldes udskillelsen af hormonet fra blodet, som bl.a. foregår i nyrerne. Udskillelsen finder faktisk sted under hele forløbet (mere om det senere), men grunden til at LH-koncentrationen først begynder at falde efter noget tid, er at produktionen af hormonet bliver mindre og mindre med tiden. Nu opstår dog spørgsmålet, hvorfor produktionen af hormonet bliver mindre og mindre med tiden. Det er lidt sværeste at forklare. Det skyldes nemlig to ting:
Jeg er dog ikke helt færdig, for hvad er mon grunden til, at LH-releasing-hormon-koncentrationen bliver mindre med tiden?
Forklaringen er, at det udskilles fra blodet gennem nyrerne hurtigere end det kan dannes.
LH-releasing-hormon dannes naturligt i Hypotalamus i hjernen, men det viser sig, at målcelle-hormonet påvirker hypotalamus til at mindske den naturlige produktion af LH-releasing-hormon (se Figur 2).
Så LH stiger, hvilket medfører mere målcelle-hormon, hvilket medfører mindre LH-releasing-hormon, hvilket medfører mindre LH osv. osv.
Det kan altså godt lade sig gøre at forstå dette system på et vist niveau uden at inddrage røre matematik. Jeg ser dog mindst tre grunde til at det alligevel er oplagt at anvende matematik på det:
Den relevante matematik er faserum, vektorfelter, og systemer af ordinære førsteordens differentialligninger. Hvis dette er ukendt for læseren, kan jeg henvise til 3b1b's fremragende video om stoffet: 3b1b: Differential equations, studying the unsolvable
Mit primære fokus i dette projekt var at opstille en matematisk model for systemet og bruge numeriske metoder til at anvende og studere modellen. Kernen i opgaven er det her faserum, som består af alle mulige tilstande af systemet. Der er tre variable, nemlig koncentrationen af hhv. LH-releasing-hormon, LH og målcellehormon, så faserummet er i 3D:
Her har jeg valgt \(s\) for "state". Der ønskes en model af, hvordan en given tilstand ændrer sig over tid for et givent system. Mere præcist ønsker vi at beskrive den hastighed - altså en 3D-vektor - en given tilstand bevæger sig med i det abstrakte faserum (3b1b forklarer rigtig godt. Se hans video!). For klarhedens skyld har jeg animeret, hvordan det evt. kunne se ud, når systemet udvikler sig og 3D-punktet - der repræsenterer systemets tilstand - flytter sig:
Den røde vektor angiver hastigheden af punktet.
Matematisk set fås denne vektor ved at differentiere tilstanden \(\vec{s}\) med hensyn til tiden:
Dette er et system af differentialligninger, men det er faktisk også et vektorfelt, for \(\large \frac{d\vec{s}}{dt}\) afhænger nemlig af \(\vec{s}\) snarere end \(t\) (Systemet har ingen hukommelse).
Grunden til at vi er så interesserede i dette vektorfelt er at det giver os mulighed for at simulere systemets udvikling over tid. Alternativt kunne man også prøve at løse differentialligningssystemet analytisk, men det er som regel umuligt.
Vi skal altså have opstillet tre differentialligniner - én pr. hormon, dvs. vi skal have fundet på nogle fornuftige modeller for, hvordan hormon-koncentrationen ændrer sig pr. tidsenhed som funktion af tilstanden \(\vec{s}\). Vi tager dem én ad gangen:
\( \large \frac{d}{dt} [LHreleasinghormon] \)
Ændringen i koncentrationen af LH-releasing-hormonet (hypothalamus-hormonet) - afhænger af hvor hurtigt hypothalamus arbejder fratrukket hvor hurtigt kroppen udskiller hormonet fra blodet (det sker vist i nyrerne). Hvor hurtigt hypothalamus arbejder er betinget af koncentrationen af hormonet dannet i målcellerne. Jo højere koncentration, jo mindre aktiv er hypothalamus. Jeg søger derfor en funktion, som er størst, når koncentrationen af målcelle-hormonet er lav og nærmer sig 0, når koncentrationen af målcelle-hormonet er høj. Eksponentialfunktionen virker som et oplagt valg:
Jeg antager desuden, at hvor hurtigt de tre hormoner udskilles fra blodet er ligefrem proportionalt med koncentrationen (Fra matematikken ved vi, at det giver en eksponentiel (negativ) vækst). En fin model for ændringen i koncentrationen af LH-releasing-hormonet (pr. tidsenhed) er altså:
Her er \(c_1\), \(c_2\) og \(c_3\) konstanter, hvor \(c_1\) kan tænkes som den maksimale produktionshastighed af LH-releasing-hormonet. \(c_2\) fortæller noget om, hvor hurtigt produktionen af LH-releasing-hormonet mindskes, når [målcellehormonet] øges. \(c_3\) bestemmer, hvor hurtigt LH-releasing-hormonet fjernes fra blodet.
\( \large \frac{d}{dt} [LH] \)
For LH-hormonet søger vi en differentialligning, som beskriver det forhold, at hypofysen reagerer med positiv feedback på LH-releasing-hormonet, og negativ feedback på målcelle-hormonet, dvs. en funktion af to variable. Det kunne fx være denne funktion:
Det ses, at produktionen af LH er lav når [målcelle-hormon] er høj. Den maksimale produktionshastighed sker når [LH-releasing-hormon] er høj, og [målcelle-hormon] er fraværende.
Funktionen er fremkommet ud fra følgende ræsonnement:
Antag først at målcelle-hormonet ikke er til stede ved hypofysen. Da er produktionen af LH alene en funktion af LH-releasing-hormonet. Det virker oplagt, at produktionen af LH er 0, når [LH-releasing-hormon] er 0. Produktionen stiger, når [LH-releasing-hormon] vokser, men den kan selvfølgelig kun vokse til et vist punkt (Der er grænser for, hvor hurtigt en kirtel kan arbejde). Jeg vil forvente en funktion, der vokser langsommere og langsommere og til sidst flader ud, og til det formål vil jeg anvende, hvad jeg vil kalde "faldskærmsfunktionen" (da man kan udlede, at det faktisk er denne v(t)-funktion der gør sig gældende ved et frit fald med luftmodstand. Det antages at luftmodstanden er kvadratisk proportional med farten). Den ser sådan her ud:
Omskrivningen viser, at det er en sum af to logistiske funktioner. Grafen for \(f\) er vist her:
Så nu har vi en model for produktionshastigheden af LH, når målcellehormonet ikke er til stede. Men hvad hvis det er til stede? Jeg fik den idé at knytte et tal mellem 0 og 1 til hver koncentation af målcellehormonet og gange det med \(f\). Det sikrer at grafens facon bevares for en fastholdt værdi af [målcellehormon]. Tallet der ganges på \(f\) skal være 1, når [målcellehormon] er 0 og nærme sig 0, når [målcellehormon] øges. Eksponentialfunktionen virker altså igen omlagt. Det giver os et rimeligt bud på modellen:
Det er egentlig sjovt, at LH-producenten har denne noget sofistikerede funktion indbygget i sig! Måske det skal være mit næste projekt at undersøge, hvordan det virker.
Her er \(c_4\) en konstant der kan fortolkes som den maksimale produktionshastighed af LH.
\(c_5\) er et tal der bestemmer hvor kraftigt LH-releasing-hormonet virker på hypofysens produktion af LH, dvs. hvor meget af stoffet der skal til for at opnå fx 50% af maksimumkapaciteten.
\(c_6\) bestemmer hvor kraftigt målcellehormonet virker på at dæmpe hypofysens produktion af LH.
Jeg har ligesom før fratrukket leddet \(c_7 \cdot [LH]\) for at tage højde for det, der fjernes af nyrerne.
\( \large \frac{d}{dt} [målcellehormon] \)
Den sidste differentialligning skal beskrive ændringen i koncentrationen af målcellehormonet pr. tidsenhed. Målcellerne reagerer på LH med positiv feedback. Jeg benytter igen "faldskærmsfunktionen" (\(f\)). Modellen bliver dermed:
Vi har altså opstillet dette vektorfelt:
Vi har nu styr på matematik-delen. Lad os nu kigge på, hvordan det kan implementeres i Python:
import numpy as np
class State:
def __init__(self, lhreleasinghormon, lh, maalcellehormon):
self.lhreleasinghormon = lhreleasinghormon
self.lh = lh
self.maalcellehormon = maalcellehormon
def get_content(self):
return [self.lhreleasinghormon, self.lh, self.maalcellehormon]
class System:
def __init__(self, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10):
self.c1 = c1
self.c2 = c2
self.c3 = c3
self.c4 = c4
self.c5 = c5
self.c6 = c6
self.c7 = c7
self.c8 = c8
self.c9 = c9
self.c10 = c10
def faldskaerm(self, x):
if x > 100:
return 1
elif x < -100:
return -1
else:
return (np.exp(x) - 1) / (np.exp(x) + 1)
def derivative(self, lhreleasinghormon, lh, maalcellehormon):
lhreleasinghormon_dot = self.c1 * np.exp(- self.c2 * maalcellehormon) - self.c3 * lhreleasinghormon
lh_dot = self.c4 * self.faldskaerm(self.c5 * lhreleasinghormon) * np.exp(- self.c6 * maalcellehormon) - self.c7 * lh
maalcellehormon_dot = self.c8 * self.faldskaerm(self.c9 * lh) - self.c10 * maalcellehormon
return lhreleasinghormon_dot, lh_dot, maalcellehormon_dot
def upcoming_state(self, state, time_step):
lhreleasinghormon_dot, lh_dot, maalcellehormon_dot = self.derivative(
state.lhreleasinghormon,
state.lh,
state.maalcellehormon
)
new_lhreleasinghormon = state.lhreleasinghormon + lhreleasinghormon_dot * time_step
new_lh = state.lh + lh_dot * time_step
new_maalcellehormon = state.maalcellehormon + maalcellehormon_dot * time_step
return State(new_lhreleasinghormon, new_lh, new_maalcellehormon)
def future_state(self, intitial_state, total_time, time_step):
current_state = intitial_state
for time in np.arange(0, total_time, time_step):
current_state = self.upcoming_state(current_state, time_step)
return current_state
def future_state_vectorized(self, intitial_state, sorted_times, time_step):
"""
Antager at sorted_times er sorteret og
består af unikke, ikke-negative tal
og at mindste tidsforskel i
sorted_times er større end time_step.
"""
current_state = intitial_state
states = []
counter = 0
time = 0
while counter < len(sorted_times):
if time >= sorted_times[counter]:
states.append(current_state)
counter += 1
current_state = self.upcoming_state(current_state, time_step)
time += time_step
return states
Så logikken i koden er basically at der er de her to klasser State og System. State er den tilstand, systemet er i, og System er selve systemet som har tilknyttet en masse funktioner, herunder future_state, som giver mulighed for at beregne fremtidige tilstande inden for systemet.
Her har jeg brugt koden til at finde frem til koncentrationen af de tre hormoner ved starttidspunktet. Jeg har lavet en genetisk algoritme, inspireret af princippet naturlig selektion, til at finde frem til de 10 konstanter i modellen ved at bruge målepunkterne fra opgaven. Her ses hvordan algoritmen langsomt kommer tættere og tættere på at beskrive dataene:
import numpy as np
import matplotlib.pyplot as plt
def gamma_ish(mu_list, sigma_list):
sample_list = []
for i in range(len(mu_list)):
sample = -1
while sample < 0:
sample = np.random.normal(mu_list[i], sigma_list[i])
sample_list.append(sample)
return np.array(sample_list)
def cost_function_LH(system, state, data_tuple, time_step):
x_data, y_data = data_tuple
container = my_system.future_state_vectorized(
init_state, x_data, time_step
)
unpacked_states = [state.get_content() for state in container]
matrix1 = np.array(unpacked_states).T
error = np.sum((y_data - matrix1[1]) ** 2)
return error, matrix1[1]
data = [
(0, 50),
(15, 338.84766),
(30, 428.71816),
(50, 400.18376),
(60, 360.62972),
(120, 208.98756),
(180, 156.81751),
]
x_data = [point[0] for point in data]
y_data = [point[1] for point in data]
system_knobs = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=float)
system_state = np.array([7, 9, 13], dtype=float)
time_step = 1
spread = 0.4
error = np.inf
output_file = "evolution_data.npz"
output_file_extra = "evolution_data_extra.npz"
fig, ax = plt.subplots()
def run_optimization():
global error, system_knobs, system_state, my_system, init_state, counter
new_system_knobs = gamma_ish(system_knobs, spread * system_knobs)
new_system_state = gamma_ish(system_state, spread * system_state)
my_system = main.System(*new_system_knobs)
init_state = main.State(*new_system_state)
new_error, evolution_extra = cost_function_LH(my_system, init_state, (x_data, y_data), time_step)
if new_error < error:
system_knobs = new_system_knobs
system_state = new_system_state
error = new_error
counter = 0
def run():
global counter
while True:
counter += 1
run_optimization()
# print(counter)
time.sleep(0.001)
run()
Algoritmen når frem til at konstanterne i systemet må være:
c1= 6.03284186e-02
c2= 5.66676837e+00
c3= 3.23887793e-03
c4= 3.38888400e+01
c5= 7.66748974e+00
c6= 9.47944174e-04
c7= 4.38393164e-02
c8= 2.20801044e+01
c9= 7.30976459e-02
c10= 1.00794645e-02