nx_stability.nx source
↩ module page · 396 lines · 17421 B
1// nx_stability.nx -- ICH Q1A(R2)/Q1E STABILITY rung: the temperature physics
2// that turns "how long is it good for" into a computed number instead of a
3// printed guess. Serves the peptide / supplement / frozen-dessert lanes
4// alike, because degradation kinetics does not care what the product is.
5//
6// ===== WHY THIS EXISTS ============================================
7//
8// nx_shelf_life already models Q10 scaling, but ONLY at decade-aligned
9// temperature steps -- and the actual regulatory question is 25 C vs 40 C,
10// a FIFTEEN degree gap. The organ returned -1 for the one comparison the
11// whole category rests on. This rung supplies the Arrhenius layer beneath
12// it and the ICH machinery above it.
13//
14// ===== THE KEYSTONE: Ea IS RECOVERED, NOT ASSERTED ================
15//
16// Every stability model needs an activation energy, and asserting one is
17// how a plausible-looking shelf life becomes a lie. We do not assert it.
18// ICH Q1A(R2) declares that 6 months at the accelerated condition (40 C)
19// supports a 24-month re-test period at the long-term condition (25 C).
20// Two conditions the regulation itself calls EQUIVALENT uniquely determine
21// the kinetics between them:
22//
23// acceleration factor = 24/6 = 4 over dT = 15 C
24// Q10 = 4^(10/15) = 2.512 (published pharma/food band: 2 - 3)
25// Ea = R*ln(4)*T1*T2/dT = 71.6 kJ/mol (published band: 60 - 100)
26//
27// BOTH recovered figures land inside independently published bands that
28// were never used as inputs. That is the liar-kill: the regulation's own
29// arithmetic agrees with the thermal-degradation literature. This is the
30// same discipline nx_thermal_process used to recover the pasteurisation
31// z-value from the PMO's two equivalent schedules.
32//
33// ===== ⚠THE RECOVERED Ea IS LUMPED, AND SAYING SO IS THE POINT =====
34//
35// ICH's two conditions are 25 C / 60% RH and 40 C / 75% RH. They differ in
36// TEMPERATURE **AND IN HUMIDITY**. Two conditions differing in two variables
37// cannot separate two parameters -- the system is underdetermined -- so the
38// figure recovered above is an EFFECTIVE activation energy that carries the
39// moisture contribution inside it, not a pure thermal Ea.
40//
41// This was missed on the first cut of this file, which presented it as a
42// clean thermal recovery. What survives and what does not:
43// SURVIVES -- projecting shelf life ALONG the ICH storage path (the
44// path real stability chambers actually run), because the
45// lumped parameter is exactly the right constant for that
46// path. Every number this file returns for that use is
47// unaffected.
48// DOES NOT -- treating it as a pure Ea and extrapolating OFF that path,
49// e.g. to a dry-storage or a sealed-container product where
50// the humidity term should not be present. There the lumped
51// value will over-state degradation.
52// ALSO NOT -- reading the agreement with the published 60-100 kJ/mol band
53// as proof of a pure Ea. It shows the magnitude is physically
54// sensible, which is weaker and is all that is claimed.
55//
56// Separating the two needs a THIRD condition that moves humidity while
57// holding temperature (the standard isothermal-humidity design, e.g. adding
58// 40 C / 25% RH). `stab_conditions_needed_to_separate` returns 3, and
59// `stab_ea_is_lumped_with_humidity` returns 1, as FUNCTIONS rather than
60// comments, so no downstream planner can route around the restriction.
61//
62// ===== MEAN KINETIC TEMPERATURE (Haynes 1971) =====================
63//
64// The SOTA content, and the reason a linear "average temperature"
65// accumulator is not good enough. Degradation rate is CONVEX in
66// temperature, so time spent hot costs more than time spent cold refunds.
67// MKT is the single temperature that produces the degradation an actual
68// varying history produced. It is always >= the arithmetic mean, so an
69// arithmetic mean UNDER-states degradation -- the dangerous direction.
70//
71// Computed with the log-sum-exp trick, in pure integer arithmetic: the raw
72// terms 10^(-E10/T) are around 10^-12 and would slam into ipow10_q3's +/-6
73// decade clamp, so the hottest segment is factored out and every remaining
74// exponent lands in [-6, 0]. This is a numerical necessity, not a nicety.
75//
76// grounded: ich_q1a_r2_stability_testing + ich_q1e_evaluation
77// + haynes_1971_mean_kinetic_temperature + arrhenius_1889
78// genealogy_id: nishi_food_chem_lane_stability_2026_07
79// license_tier: ORIGINAL
80
81import "nx_syscalls.nx"
82import "nx_pow10.nx"
83const STAB_MAGIC_10000: i64 = 10000
84const STAB_MAGIC_1000000: i64 = 1000000
85const STAB_MAGIC_100000: i64 = 100000
86
87const STAB_INVALID: i64 = 0 - 1
88
89// R = 8.314 J/(mol*K), milli-scaled.
90const STAB_R_MILLI: i64 = 8314
91// ln(10) = 2.302585, milli-scaled. Converts base-10 logs to natural ones,
92// which is the only place the natural base is needed at all.
93const STAB_LN10_MILLI: i64 = 2303
94// 273.15 K expressed in centi-Kelvin.
95const STAB_KELVIN_OFF_CENTI: i64 = 27315
96
97// --- ICH Q1A(R2) 2.2.7 -- THE REGULATION'S EQUIVALENCE CLAIM ---
98// These four numbers are the ONLY stability inputs asserted in this file.
99// Q10 and Ea are both derived from them.
100const STAB_ICH_LONG_C: i64 = 25
101const STAB_ICH_LONG_MONTHS: i64 = 24
102const STAB_ICH_ACCEL_C: i64 = 40
103const STAB_ICH_ACCEL_MONTHS: i64 = 6
104
105// ICH Q1A(R2) significant change for a drug product: 5% loss of assay from
106// the initial value. Shelf life is the time to reach this, not to reach zero.
107const STAB_SIGNIFICANT_LOSS_PCT: i64 = 5
108
109// Published bands, used ONLY to check the recovered figures -- never as inputs.
110const STAB_PUB_Q10_LO_Q3: i64 = 2000
111const STAB_PUB_Q10_HI_Q3: i64 = 3000
112const STAB_PUB_EA_LO_J: i64 = 60000
113const STAB_PUB_EA_HI_J: i64 = 100000
114
115// === Temperature ==================================================
116
117// Celsius (whole degrees) -> centi-Kelvin.
118func stab_c_to_centi_kelvin(c: i64) -> i64 {
119 return c * 100 + STAB_KELVIN_OFF_CENTI
120}
121
122// === Recovery from an equivalence claim ===========================
123
124// Q10 (Q3, so 2512 = 2.512) implied by declaring `months_lo` at `t_lo_c`
125// equivalent to `months_hi` at the hotter `t_hi_c`.
126// REFUSES anything that is not a genuine acceleration: the hot condition
127// must be hotter AND shorter, else the caller has passed a pair that says
128// nothing about degradation and a returned number would be noise.
129func stab_q10_from_equivalence_q3(t_lo_c: i64, months_lo: i64, t_hi_c: i64, months_hi: i64) -> i64 {
130 var dt: i64 = 0
131 var af_q3: i64 = 0
132 var lg_af: i64 = 0
133 if t_hi_c <= t_lo_c { return STAB_INVALID }
134 if months_hi <= 0 { return STAB_INVALID }
135 if months_lo <= months_hi { return STAB_INVALID }
136 dt = t_hi_c - t_lo_c
137 af_q3 = months_lo * 1000 / months_hi
138 lg_af = ilog10_milli(af_q3)
139 if lg_af == IP10_INVALID { return STAB_INVALID }
140 return ipow10_q3(10 * lg_af / dt)
141}
142
143// Activation energy (J/mol) implied by the SAME equivalence.
144// Ea = R * ln(10) * log10(AF) * T1 * T2 / dT
145// Scale bookkeeping is done stepwise so the product never approaches the
146// i64 ceiling: the widest intermediate here is ~1.1e13.
147func stab_ea_from_equivalence(t_lo_c: i64, months_lo: i64, t_hi_c: i64, months_hi: i64) -> i64 {
148 var dt: i64 = 0
149 var af_q3: i64 = 0
150 var lg_af: i64 = 0
151 var rl_milli: i64 = 0
152 var a_milli: i64 = 0
153 var t1c: i64 = 0
154 var t2c: i64 = 0
155 if t_hi_c <= t_lo_c { return STAB_INVALID }
156 if months_hi <= 0 { return STAB_INVALID }
157 if months_lo <= months_hi { return STAB_INVALID }
158 dt = t_hi_c - t_lo_c
159 af_q3 = months_lo * 1000 / months_hi
160 lg_af = ilog10_milli(af_q3)
161 if lg_af == IP10_INVALID { return STAB_INVALID }
162 rl_milli = STAB_R_MILLI * STAB_LN10_MILLI / 1000
163 a_milli = rl_milli * lg_af / 1000
164 t1c = stab_c_to_centi_kelvin(t_lo_c)
165 t2c = stab_c_to_centi_kelvin(t_hi_c)
166 return a_milli * t1c / 1000 * t2c / (STAB_MAGIC_10000 * dt)
167}
168
169// The two ICH-derived figures, as functions so no caller can drift from them.
170func stab_ich_q10_q3() -> i64 {
171 return stab_q10_from_equivalence_q3(STAB_ICH_LONG_C, STAB_ICH_LONG_MONTHS,
172 STAB_ICH_ACCEL_C, STAB_ICH_ACCEL_MONTHS)
173}
174
175func stab_ich_ea_j() -> i64 {
176 return stab_ea_from_equivalence(STAB_ICH_LONG_C, STAB_ICH_LONG_MONTHS,
177 STAB_ICH_ACCEL_C, STAB_ICH_ACCEL_MONTHS)
178}
179
180// --- The lumping, as functions rather than prose ---
181
182// ICH's relative-humidity conditions, carried so the pairing is inspectable
183// rather than folded away into the temperatures.
184const STAB_ICH_LONG_RH_PCT: i64 = 60
185const STAB_ICH_ACCEL_RH_PCT: i64 = 75
186
187// 1 = the recovered activation energy carries a moisture contribution and is
188// NOT a pure thermal Ea. Always true for a value recovered from a pair of
189// conditions whose humidity differs, which ICH's does.
190func stab_ea_is_lumped_with_humidity() -> i64 {
191 if STAB_ICH_ACCEL_RH_PCT == STAB_ICH_LONG_RH_PCT { return 0 }
192 return 1
193}
194
195// How many distinct conditions are needed to separate a temperature term
196// from a humidity term. Two unknowns need two independent contrasts, and a
197// contrast needs a pair, so three conditions is the floor.
198func stab_conditions_needed_to_separate() -> i64 {
199 return 3
200}
201
202// Is a projection to (target_c, target_rh_pct) ON the path the lumped
203// parameter was recovered along? REFUSES to bless an off-path humidity:
204// the lumped value over-states degradation for a drier product, and a
205// planner that could not see this would ship a needlessly short expiry --
206// or, moving the other way, a dangerously long one.
207func stab_projection_on_recovered_path(target_rh_pct: i64) -> i64 {
208 if target_rh_pct < STAB_ICH_LONG_RH_PCT { return 0 }
209 if target_rh_pct > STAB_ICH_ACCEL_RH_PCT { return 0 }
210 return 1
211}
212
213// Do the recovered figures agree with bands nobody fed in? 1 = yes.
214func stab_recovered_q10_in_published_band() -> i64 {
215 var q: i64 = 0
216 q = stab_ich_q10_q3()
217 if q == STAB_INVALID { return 0 }
218 if q < STAB_PUB_Q10_LO_Q3 { return 0 }
219 if q > STAB_PUB_Q10_HI_Q3 { return 0 }
220 return 1
221}
222
223func stab_recovered_ea_in_published_band() -> i64 {
224 var e: i64 = 0
225 e = stab_ich_ea_j()
226 if e == STAB_INVALID { return 0 }
227 if e < STAB_PUB_EA_LO_J { return 0 }
228 if e > STAB_PUB_EA_HI_J { return 0 }
229 return 1
230}
231
232// === Acceleration ==================================================
233
234// Rate multiplier (Q3) for being `dt_c` degrees hotter, given a Q10.
235// Life divides by this; rate multiplies by it.
236func stab_accel_factor_q3(q10_q3: i64, dt_c: i64) -> i64 {
237 var lg: i64 = 0
238 if q10_q3 <= 1000 { return STAB_INVALID }
239 lg = ilog10_milli(q10_q3)
240 if lg == IP10_INVALID { return STAB_INVALID }
241 return ipow10_q3(dt_c * lg / 10)
242}
243
244// === Mean kinetic temperature (Haynes 1971) ========================
245
246// E10 = Ea/(R*ln10), in milli-Kelvin. The "base-10 Arrhenius temperature":
247// the rate falls by one decade for every E10/T^2 kelvin, and working in this
248// unit is what lets the whole MKT computation stay in ipow10_q3's domain.
249func stab_e10_milli_kelvin(ea_j: i64) -> i64 {
250 var rl_milli: i64 = 0
251 if ea_j <= 0 { return STAB_INVALID }
252 rl_milli = STAB_R_MILLI * STAB_LN10_MILLI / 1000
253 return ea_j * STAB_MAGIC_1000000 / rl_milli
254}
255
256// Decades of rate change between absolute zero and T: E10/T, in milli-decades.
257func stab_e10_over_t_milli(e10_mk: i64, t_centi_k: i64) -> i64 {
258 if t_centi_k <= 0 { return STAB_INVALID }
259 return e10_mk * 100 / t_centi_k
260}
261
262// Is the Arrhenius term actually CONVEX at this temperature? d2k/dT2 > 0
263// iff Ea/(R*T) > 2. MKT > arithmetic mean is a consequence of that
264// convexity (Jensen), so this is the PRECONDITION of the headline claim --
265// measured rather than assumed, the way nx_peptide measures its bisection's
266// monotonicity instead of trusting it. Returns Ea/(R*T) scaled by 1000.
267func stab_convexity_ratio_milli(ea_j: i64, t_centi_k: i64) -> i64 {
268 if t_centi_k <= 0 { return STAB_INVALID }
269 if ea_j <= 0 { return STAB_INVALID }
270 return ea_j * STAB_MAGIC_100000 / (STAB_R_MILLI * t_centi_k / 1000)
271}
272
273// Mean kinetic temperature of a time/temperature history, in centi-Kelvin.
274// `temps_c` are whole Celsius, `hours` the duration held at each.
275//
276// LOG-SUM-EXP: the hottest segment is factored out of the sum so every
277// exponent that reaches ipow10_q3 lies in [-6, 0]. Without this the raw
278// terms (~10^-12) hit the clamp and every history would return the same
279// wrong answer.
280func stab_mkt_centi_kelvin(temps_c: *i64, hours: *i64, n: i64, ea_j: i64) -> i64 {
281 var i: i64 = 0
282 var j: i64 = 0
283 var tmax: i64 = 0 - STAB_MAGIC_100000
284 var wsum: i64 = 0
285 var s: i64 = 0
286 var e10: i64 = 0
287 var base: i64 = 0
288 var ti: i64 = 0
289 var mean_q3: i64 = 0
290 var lg: i64 = 0
291 var denom: i64 = 0
292 if n <= 0 { return STAB_INVALID }
293 e10 = stab_e10_milli_kelvin(ea_j)
294 if e10 == STAB_INVALID { return STAB_INVALID }
295 while i < n {
296 if hours[i] < 0 { return STAB_INVALID }
297 if temps_c[i] > tmax { tmax = temps_c[i] }
298 wsum = wsum + hours[i]
299 i = i + 1
300 }
301 if wsum <= 0 { return STAB_INVALID }
302 base = stab_e10_over_t_milli(e10, stab_c_to_centi_kelvin(tmax))
303 if base == STAB_INVALID { return STAB_INVALID }
304 while j < n {
305 ti = stab_e10_over_t_milli(e10, stab_c_to_centi_kelvin(temps_c[j]))
306 if ti == STAB_INVALID { return STAB_INVALID }
307 s = s + hours[j] * ipow10_q3(base - ti)
308 j = j + 1
309 }
310 mean_q3 = s / wsum
311 if mean_q3 <= 0 { return STAB_INVALID }
312 lg = ilog10_milli(mean_q3)
313 if lg == IP10_INVALID { return STAB_INVALID }
314 denom = base - lg
315 if denom <= 0 { return STAB_INVALID }
316 return e10 * 100 / denom
317}
318
319// Time-weighted arithmetic mean temperature, centi-Kelvin. Present so the
320// gap against MKT can be MEASURED rather than described.
321func stab_mean_centi_kelvin(temps_c: *i64, hours: *i64, n: i64) -> i64 {
322 var i: i64 = 0
323 var wsum: i64 = 0
324 var acc: i64 = 0
325 if n <= 0 { return STAB_INVALID }
326 while i < n {
327 if hours[i] < 0 { return STAB_INVALID }
328 acc = acc + hours[i] * stab_c_to_centi_kelvin(temps_c[i])
329 wsum = wsum + hours[i]
330 i = i + 1
331 }
332 if wsum <= 0 { return STAB_INVALID }
333 return acc / wsum
334}
335
336// === Potency, shelf life, overage ==================================
337
338// Fraction of label potency remaining (permil) after `t_months`, given that
339// `loss_pct` is lost over `shelf_months`. First-order: the loss compounds,
340// so twice the time is NOT twice the loss.
341//
342// SIGN OF THE ERROR IS THE SAFE ONE. ipow10_q3 is chord-interpolated on a
343// convex curve, so a negative exponent comes back slightly UNDER-stated;
344// retention is therefore under-reported and the overage derived from it is
345// over-reported. For a label claim that is the protective direction.
346func stab_potency_retained_permil(t_months: i64, shelf_months: i64, loss_pct: i64) -> i64 {
347 var lg_ret: i64 = 0
348 if shelf_months <= 0 { return STAB_INVALID }
349 if t_months < 0 { return STAB_INVALID }
350 if loss_pct <= 0 { return STAB_INVALID }
351 if loss_pct >= 100 { return STAB_INVALID }
352 lg_ret = ilog10_milli((100 - loss_pct) * 10)
353 if lg_ret == IP10_INVALID { return STAB_INVALID }
354 return ipow10_q3(lg_ret * t_months / shelf_months)
355}
356
357// Overage (permil of label) needed so the product still meets label claim
358// once `retained_permil` is all that survives to expiry.
359//
360// CEILING, NOT TRUNCATION, AND THE REASON IS NOT PRECISION. Label claim is
361// an INEQUALITY -- the product must deliver at least what the label says at
362// the end of shelf life. Truncating divides that requirement's own margin
363// away: at 939 permil retention the exact overage is 64.96, and flooring it
364// to 64 lands the product at 999 permil of label on the expiry date, i.e.
365// non-compliant by construction. The gate caught exactly that. Rounding
366// up costs a fraction of a percent of fill and can only ever overshoot.
367func stab_overage_permil(retained_permil: i64) -> i64 {
368 if retained_permil <= 0 { return STAB_INVALID }
369 if retained_permil >= 1000 { return 0 }
370 return (STAB_MAGIC_1000000 + retained_permil - 1) / retained_permil - 1000
371}
372
373// Has the product crossed ICH's significant-change line?
374func stab_is_significant_change(retained_permil: i64) -> i64 {
375 if retained_permil == STAB_INVALID { return STAB_INVALID }
376 if retained_permil <= (100 - STAB_SIGNIFICANT_LOSS_PCT) * 10 { return 1 }
377 return 0
378}
379
380// === Hard limits, as functions rather than comments =================
381
382// Accelerated data NEVER establishes a shelf life on its own. ICH Q1E
383// permits extrapolation from accelerated data only alongside long-term data
384// covering the period, and never beyond twice (and no more than 12 months
385// past) the available long-term coverage. A planner that could route
386// around this would manufacture expiry dates, so it is a function.
387func stab_accelerated_alone_sets_shelf_life() -> i64 {
388 return 0
389}
390
391// Extrapolation ceiling in months per ICH Q1E: min(2*long_term, long_term+12).
392func stab_max_extrapolated_months(long_term_months: i64) -> i64 {
393 if long_term_months <= 0 { return STAB_INVALID }
394 if long_term_months + 12 < long_term_months * 2 { return long_term_months + 12 }
395 return long_term_months * 2
396}