Maternal Transfer and Long-Term Population Effects of PCBs in Baltic Grey Seals Using a New Toxicokinetic–Toxicodynamic Population Model

Empirical evidence has shown that historical exposure of polychlorinated biphenyls (PCBs) to Baltic grey seals not only severely affected individual fitness, but also population growth rates and most likely caused the retarded recovery rate of the depleted population for decades. We constructed a new model which we term a toxicokinetic–toxicodynamic (TKTD) population model to quantify these effects. The toxicokinetic sub-model describes in detail the bioaccumulation, elimination and vertical transfer from mother to offspring of PCBs and is linked to a toxicodynamic model for estimation of PCB-related damage, hazard and stress impacts on fertility and survival rates. Both sub-models were linked to a Leslie matrix population model to calculate changes in population growth rate and age structure, given different rates of PCB exposure. Toxicodynamic model parameters related to reproductive organ lesions were calibrated using published historical data on observed pregnancy rates in Baltic grey seal females. Compared to empirical data, the TKTD population model described well the age-specific bioaccumulation pattern of PCBs in Baltic grey seals, and thus, the toxicokinetic parameters, deduced from the literature, are believed to be reliable. The model also captured well the general effects of PCBs on historical population growth rates. The model showed that reduced fertility due to increased PCB exposure causes decreased vertical transfer from mother to offspring and in turn increased biomagnification in non-breeding females. The developed TKTD model can be used to perform population viability analyses of Baltic grey seals with multiple stressors, also including by-catches and different hunting regimes. The model can also be extended to other marine mammals and other contaminants by adjustments of model parameters and thus provides a test bed in silico for new substances. Supplementary Information The online version contains supplementary material available at 10.1007/s00244-022-00962-3.


Model Overview and Calculation Scheme
Our matrix model, called the TKTD population model, analyses adverse effects of PCBs on Baltic grey seals by linking toxicokinetics (TK) of dietary and vertically transferred PCB accumulation to toxicodynamics (TD) of adverse effects on reproduction and survival that ultimately predict changes to population demographics and growth rates. Survival and fecundity parameters of an age-structured Leslie matrix model for an ideal Baltic grey seal population are modified due to internal PCB concentrations following DEBtox and TDM principles of damage, hazard, stress, and recovery.
While the model runs on an annual basis for the population dynamics, each year is divided into three periods (lactation, delay and gestation), representing different lifehistory phases for a reproducing grey seal female (Fig. 1). The year starts with the lactation period (lasting for 18 days), during which pups are nursed by their mothers. After weaning, females mate and the delay period starts (lasting for 100 days), after which the fertilized egg is implanted. During the gestation period (lasting for 247 days), the fetus develops in the uterus. The year ends with the birth of new pups. Since all deaths are assumed to occur just before the birth events, all pups that are born have mothers that nurse them. Different TD effects of PCBs occur after specific time periods (e.g., some fertility effects are estimated after the delay period). The TKTD population model consist of three interconnected sub models (Fig. 2). The toxicokinetic (TK) model of PCB accumulation follows the total body PCB concentrations in seals of different age classes at the end of each time period of a year. The toxicodynamic (TD) model follows damage, hazard and reproductive stress in seals of different age classes at the end of each time period of a year. The Leslie matrix model follows number of individuals, survival rate and fertility in each age class at the end of each year. The TKTD population model was implemented as a set of function and script files in the numerical software MATLAB ® (MathWorks Inc., Natick, MA, USA) and all analyses were performed with these.
Before any calculations can be performed, a number of initial conditions have to be specified. These are initial PCB concentrations, initial damage levels and initial population vector (Table 1). Also, PCB concentrations in different prey are specified for each year of the simulation (Fig. 8). Thereafter, PCB concentrations and damage levels in fetuses, pups and females at the end of the different time periods within a year are successively calculated. Hazard/stress levels are calculated based on obtained damage levels and fertilities and survival rates are updated based on these. A new population vector is finally calculated, including the number of individuals in different age classes at the end of the year. Outputs from calculations for one year are used as inputs to calculations for the next year. The procedure is repeated for one year at a time and may be continued for as many years as desired (Fig. 2).

Fig. 2
Overview of the TKTD population model. Text boxes represent sub-models, performing calculations for year + 1. Free texts are inputs and outputs, with arrows showing relations to sub-models.
For each year, the following calculations are performed: 1. Calculation of PCB levels in pups and females at the end of the lactation period, as a result of vertical transfer from females to pups, based on PCB levels in fetuses, pups and females at the end of previous year. Transfer levels are dependent on fertilities, determining the mean number of pups per female.
2. Calculation of cumulated damage levels in pups and females at the end of the lactation period, based on PCB and damage levels in fetuses, pups and females at the end of previous year.
3. Calculation of accumulated PCB levels in pups and females at the end of the delay period, based on PCB levels at the end of the lactation period. PCB is assimilated through dietary uptake and eliminated through metabolic transformation and fecal egestion.
4. Calculation of cumulated damage levels in pups and females at the end of the delay period, based on PCB and damage levels in pups and females at the end of the lactation period.
5. Calculation of reproductive stress levels in females at the end of the delay period, based on reproductive damage levels at the end of the lactation period.
6. Calculation of reduced fertilities at the end of the delay period, based on ideal fertilities and female reproductive stress levels at the end of the delay period.
7. Calculation of PCB levels in fetuses, pups and females at the end of the gestation period, as a result of bioaccumulation and placental transfer from females to fetuses, based on PCB levels in pups and females at the end of the delay period. Transfer levels are dependent on fertilities at the end of the delay period, determining the mean number of produced fetuses per female.
8. Calculation of cumulated damage levels in fetuses, pups and females at the end of the gestation period, based on PCB and damage levels in pups and females at the end of the delay period.
9. Calculation of hazard/stress levels in fetuses, pups and females at the end of the gestation period, based on damage levels in pups and females at the end of the delay period.
10. Calculation of reduced survival rates at the end of the year, based on ideal survival rates and hazard levels in pups and females at the end of the gestation period.
11. Calculation of reduced fertilities at the end of the year, based on fertilities at the end of the delay period (reduced due to reproductive stress) and fetal hazard levels at the end of the gestation period.
12. Calculation of new population vector (accounting for population growth during one year), based on updated fertilities and survival rates, and previous year's population vector.
In order to describe the temporal dynamics of variables for fetuses, pups and females of different age classes at different times within a year, we use the notation , ( ). The first letter represents the variable of interest, here being seal PCB concentration in mg/kg. The upper index indicates the time period (l = lactation, d = delay, g = gestation). The first lower index indicates if it is a fetus ( = 0), a pup ( = 1) or a female ( = ), where = 2, … , 46 indicates the age class. Individuals of age class = 1 (with an age of 0-1 year) are here referred to as pups, whereas individuals of age class = 2, … , 46 are referred to as females. For fetuses and pups, the second lower index indicates the age class of the mother. Hence, 0, ( ), 1, ( ) and , ( ) are PCB concentrations at time during time period for fetuses with mothers of age class , pups with mothers of age class and females of age class , respectively. The temporal dynamics is described by differential equations. By solving these, state variables at the end of time periods within a year are calculated. These are used to perform step-wise calculations from one time period to the next and also from one year to the next, using obtained variable values as initial conditions for the next calculation. In order to capture variable states at the end of different time periods for specified years , we also introduce the notations ̂0 ,i ( ), ̂1 ,i ( ) and ̂( ), where a hat over a variable name indicates reference to the end of a time period. For example, ̂1 , ( ) are PCB concentrations at the end of the lactation period (l) of year in pups with mothers of age class . The same conventions are used also for other variables (damages, hazards and stresses). Furthermore, we use the notation where a dot above letters indicates a rate with respect to time and a straight overline indicates an average value. The state variables of the TKTD model are listed in Table 1.

Toxicokinetic Model
The toxicokinetic model describes PCB bioaccumulation in pups and females through dietary uptake and elimination as well as vertical transfer from mother to fetus through the placenta during gestation and from mother to pup through breast milk during lactation (Fig. 3).

Fig. 3
The toxicokinetic model and its connections to other sub models. Text boxes represent sub modules. Free texts are inputs and outputs, with arrows showing relations to sub modules. The flowchart illustrates computations performed for year + 1, based on inputs from previous year .
During the lactation period, PCB is transferred from females to pups through the breast milk. Neither females nor pups consume prey during this period. Hence, no PCB is accumulated from the environment, but it is redistributed from females to pups through vertical transfer, including some losses to the environment. PCB concentrations also change due to growth dilution, especially in the fast growing pups. During the delay period, both females and pups feed on prey. Hence, PCB concentrations change due to dietary uptake, elimination (fecal egestion and metabolic transformation) and growth dilution. During the gestation period, females and pups continue to feed, but females also transfer PCB to fetuses through the placenta. Hence, concentration changes in pups are governed by dietary uptake, elimination, and growth dilution, whereas concentration changes in females are governed by dietary uptake, elimination, growth dilution and vertical transfer. Furthermore, the PCB concentration increases in fetuses due to vertical transfer, which affects the PCB concentration in new-born pups at the start of the next year.

Growth model
In order to calculate growth dilution rates and keep track of body weights (affecting rates of feeding, vertical transfer and PCB elimination), a sub model for body growth was established (Fig. 4).
Fetuses are assumed to growth exponentially during the gestation period (Peters 1983). Hence, the body weight of a fetus at time during a year is obtained as: Here, 00 is the initial fetal weight (the weight of a fertilized egg) and 0 is an exponential growth rate constant. At the end of the gestation period, the fetus has reached the birth weight ℎ . Hence, the initial fetal weight 00 is calculated as: For simplicity it is assumed that new-born pups continue to grow exponentially with the same growth rate as fetuses during the lactation period. The body weight of a pup at time during the lactation period is obtained as: At the end of the lactation period, the pup has reached the weaning weight; 1 (∆ ) = . Hence, the exponential growth rate constant 0 is calculated as: After weaning, grey seal pups typically lose weight during some months due to reduced nutritional intake (Boyd 1984), but when they become more skilled hunters they start to regain weight. For simplicity, it is here assumed that pups keep the weaning weight until they reach sub adult age = 1 yr.
Grey seal females are typically full-grown at an age of six years. When they nurse pups (during the lactation period) they lose a considerable amount of weight, but this is regained when they feed during the remaining part of the year. Under constant environmental conditions, dynamic energy budget (DEB) theory predicts body growth according to the von Bertalanffy growth function (Kooijman 2000), in which the body weight ( ) at age is obtained from the equation: Here, ∞ is the asymptotic body weight, 0 is the initial body weight (the body weight at = 0 ) and is a growth rate constant. The von Bertalanffy growth function is here adopted to describe body growth of grey seal females from weaning ( = ∆ ) up to the age of sexual maturation ( = 5 yr). At that age, females reach the mature body weight and start to breed. The body weight of non-mature females ( = 2 − 5) at time during the growth period of a year is obtained as: The weaning body weight , the mature body weight , the post-lactation body weights ( ,6 and ,7 ) and the maximum body weight are specified in Table 2, based on data for full-grown Baltic grey seal females (Härkönen 2016). Females lose significant weight during lactation and it is assumed that their weight decrease exponentially throughout the lactation period, from the initial body weight ( or ) to the post-lactation body weight; ,6 or ,7 for primiparous females (first-time breeders) and multiparous females (experienced breeders), respectively: Here, is the exponential body weight decline rate constant, calculated as: After lactation, during the delay period, mature females are assumed to linearly regain the weight they lost during lactation: 6 ( ) = ,6 + 6 ( − ∆ ) ( ) = ,7 + 7 ( − ∆ ) ( = 7, … , 46) The regain rate constants ( 6 and 7 ) are calculated as: After the delay period, mature females continue to grow according to von Bertalanffy throughout the gestation period and reach the maximum body weight at the end of the year: In succeeding years, females are assumed to repeat the pattern of exponential decline (lactation period), linear regain (delay period) and von Bertalanffy growth (gestation period), see Fig. 4.
The mean body weight ̅ during time period is calculated as the integrated mean value over the period ( 0 ≤ ≤ 0 + ∆ ):

Bioaccumulation model
The major routes for uptake, transfer and elimination of PCB in females, pups and fetuses are illustrated in Fig. 5. The corresponding toxicokinetic equations are now treated for each groups separately.

Bioaccumulation in females
Bioaccumulation in females older than one year ( = 2, … , 46) includes dietary uptake, elimination and vertical transfer to offspring (Fig. 5). The toxicokinetics is described by a physiologically based equation similar to the one-compartment model, commonly used in DEBtox models, but extended to account also for growth dilution and vertical transfer. The rate of change of the total amount of PCB in a female is the rate by which PCB is assimilated (assumed to be proportional to the concentration in the prey) subtracted by the rate by which PCB is eliminated (assumed to be proportional to the concentration in the body): Here, , ( ) is the total body PCB concentration (mg/kg) at time (yr) in a female of age class during time period , , is the mean PCB concentration in the diet (mg/kg) and ( ) is the female body weight (kg). The repeated notation of index in female concentrations ( , ) is applied for technical reasons, enabling a compact general description of all age classes and fetuses. The dietary uptake rates (yr -1 ) for females (and pups) are defined as: Here, ̇, is the prey consumption rate (kg/yr). The fraction of consumed PCB that is assimilated by the seal is assumed to equal the diet assimilation efficiency , the fraction of ingested prey that is assimilated into energy. The assumption is motivated by studies indicating that the assimilation of hydrophobic substances is closely linked to digestion of associated food (Bodiguel et al. 2009). Dietary uptake rates are zero when = since neither females nor pups feed during lactation. As a simplification, it is assumed that also non-mature females ( = 2, … , 5) fast during the lactation period, though they have no pups to nurse. The simplification will only have minor effects since the lactation period is very short (18 days) compared to a whole year. All seals are assumed to feed throughout the whole delay period, neglecting the two weeks of fasting during molting. According to DEB theory, the feeding rate of an animal in a constant environment is proportional to the surface area of the structural body volume (Kooijman 2000). By assuming this, isometric growth and body weight proportional to body volume, the prey consumption rate ̇ for a seal of weight is expressed as: The prey consumption factor is a constant, estimated from empirical data of prey consumption in Baltic grey seals.
Mature females lactate and hence they lose and regain weight during a year. Due to their rapid weight regain after lactation, they can be expected to peek their daily prey consumption during the delay period, though they have higher weight during other periods. In a DEB context, the reserve (blubber) mass is reduced during lactation, but the structural body mass is not. The weight in Eq. (15) is best considered as a mean body weight for females of a certain body length (which do not reduce during fasting). As a simplification, it is here assumed that the prey consumption rate for seals of age class can be described by Eq. (15) with as the maximum body weight of a year, i.e. the body weight at the end of the year ( = 1 yr). The age-specific dietary uptake rate constants for females (and pups) are approximated as: Here, ̅ is the mean body weight of individuals in age class during time period .
The female PCB elimination rates (yr -1 ) accounting for all elimination routes, are defined as: Elimination through vertical transfer is accounted for via the lactational transfer rate constant , and the placental transfer rate constant , , respectively. There are two separate routes for PCB removal through fecal egestion. The first route is accounted for via the diet assimilation efficiency ( ), with elimination proportional to the concentration in the prey. The second route is accounted for via the fecal elimination rate constant ( , ), with elimination proportional to the concentration in the body. However, the contribution from the second route is probably minor compared to the first one (Bodiguel et al. 2009). For simplicity, the removal rate constant , accounts for the combined effect of metabolic transformation and fecal elimination: Here, , is the metabolic transformation rate constant and , is the fecal elimination rate constant. The vertical transfer rate constants ( , and , ) are considered as mean values for each period. This simplification neglects that vertical transfer rates increase with the weight of the offspring (fetus or pup), which grows and assimilates increasing amounts of nutrients, resulting in vertical transfer rates that increase with time.
Notice that all rate constants are considered as be age-specific, though they may have same or similar values between age classes. For the sake of generality, Eq. (13) accounts for all assimilation and elimination routes, including vertical transfer from female to offspring, though all these will not occur at the same time. When applying Eq. (13) to different time periods of a year, different rate constants appear according to Eqs. (16) and (17).
The toxicokinetic equation (Eq. (13)) can be rewritten as the following differential equation for the PCB concentration , ( ) at time in a female of age class during time period : The growth dilution rate ( , ), the relative increase in body weight per unit time, was here introduced: It is assumed that the mean PCB concentration in the prey is constant throughout each year; , ( ) = ̅ , ( ), where ̅ , ( ) is the mean PCB concentration in the diet of seals in age class during year . Moreover, it is assumed that the growth dilution rate is approximately constant within each period of a year. The growth dilution rate ( , ) for seals of age class during a time period ( 0 ≤ ≤ 0 + ∆ ) is approximated as the mean growth dilution rate during that period: During the lactation period, body weights of pups and females are described by exponential functions (Eqs. (3) and (7)). The growth dilutions rates then simplify into the exponential growth rate constants: During the rest of the year, pups have constant body weight and hence no growth dilution: The effective elimination rate constant for females (and pups) during a time period within a year is defined as: Now, the toxicokinetic equation (Eq. (19)) can be approximated as the following first-order ordinary differential equation (ODE): The solution to Eq. (25) is: Equation (26) expresses the PCB concentration in a female of age class at time during period within a year. The first term accounts for dietary uptake (proportional to dietary PCB concentration, characterized by constant , ) and the second term accounts for elimination (exponential decrease, characterized by elimination rate ̃, ). Equation (26) is used to obtain PCB concentrations in females at the end of different time periods , based on initial PCB concentrations , (0) at the start of the periods ( = 0). Within a year, the concentration obtained in an age class at the end of a period, is used as initial value the next period. Between years, the concentration obtained in an age class ( ) at the end of a year ( ), is used as initial value in next age class ( + 1) at the start of next year ( + 1).
The PCB concentration in females of age class + 1 at the end of the lactation period ( = , = ∆ ) during year + 1 is obtained from Eq.(26) as: Here, ̂( ) is the PCB concentration in females of age class at the end of year . The mean elimination rate (accounting for lactational transfer and growth dilution) during the lactation period year is: The mean lactational transfer rate during the lactation period year + 1 is calculated as: Here, ( ) is the fertility of females in age class at the end of year . Notice that the factor 2 ( ) is the mean number of pups that each female (of age class + 1) gave birth to at the end of previous year ( ). The factor is included to account for absence of lactation in females without pups.
The PCB concentration in females of age class + 1 at the end of the delay period ( = , = ∆ ) during year + 1 is obtained from Eq. (26) as: The PCB concentration in females of age class + 1 at the end of the gestation period ( = , = ∆ ) year + 1 is obtained from Eq. (26) as: The mean elimination rate (accounting for removal, placental transfer and growth dilution) during the gestation period is: The mean placental transfer rate (for females of age class + 1) during the gestation period (year + 1) is: Here, ( ) is the fertility for females of age class at the end of the delay period year , reduced from ideal fertility due to reproductive stress (obtained from the toxicodynamic model). The factor 2 ( ) is an estimation of the mean number of fetuses that each female (of age class ) produces during year . The factor is included to account for absence of placental transfer in non-pregnant females. Equations (27), (30) and (31) are used to successively calculate PCB concentrations in females at the end of the lactation period, the delay period and the gestation period for one year at a time.

Bioaccumulation in Fetuses
Fetuses assimilate PCB only through vertical transfer via the placenta. It is assumed that all PCB that is eliminated by a mother through placental transfer is assimilated by her fetus, which has no capacity to eliminate it. The toxicokinetic equation is a mass balance where the rate of change of PCB in the fetus equals the rate at which PCB is eliminated from the mother through placental transfer. The rate of change of the total amount of PCB in a fetus with a mother of age class during the gestation period ( = ) is: Here, 0, and , are PCB concentrations in fetus and mother, 0 and are body weights of fetus and mother, whereas , is the placental transfer rate constant. The toxicokinetic equation (Eq. (34)) can be rewritten as: The fetus grows exponentially according to Eq. (1). Hence, the growth dilution rate ,0 equals the exponential growth rate: The body weight ratio between female and fetus ( 0, ) is approximated as a mean value over the gestation period: The PCB concentration in the mother ( , ) during the gestation period ( = ) is obtained from Eq. (26) and inserted into the toxicokinetic equation (Eq. (35)), which yields the following first order ODE: Equation (39) expresses the PCB concentration in a fetus with a mother of age class at time during the gestation period ( = ). The first term accounts for dietary uptake in the mother, transferred to the fetus. The second term accounts for reduction of transferred PCB due to elimination in the mother (exponential decrease, characterized by elimination rate ̃, ). The last term accounts for growth dilution in the fetus (exponential decrease, characterized by growth dilution rate ,0 ). From Eq. (39), with = ∆ , the PCB concentration in fetuses with mothers of age class at the end of the gestation period year + 1 is obtained as: Equation (40) is used to calculate PCB concentrations in fetuses at year ends, based on known PCB levels in their mothers at the end of the delay period.

Bioaccumulation in Pups
Bioaccumulation of PCB in pups (age class = 1) is described by a toxicokinetic equation, similar to Eq. (13) for females, but with no placental transfer and a positive contribution from vertical transfer through lactation. The rate of change of the total amount of PCB in a pup with a mother of age class is: Here, 1, and , are PCB concentrations in pups and mothers, ,1 is the PCB concentration in the diet, whereas 1 and are body weights of pups and mothers. The lactational transfer rate 1, for pups with mothers of age class , accounting for PCB uptake through breast milk during lactation, is: Here, , is the lactational transfer rate constant for females of age class . The fraction of transferred PCB that is assimilated by the pup is assumed to equal the lactation assimilation efficiency , the fraction of consumed milk that is assimilated into energy by a pup during lactation. The dietary uptake rate for pups ,1 is calculated analogously to females (Eq. (14)), with prey consumption rate ̇, 1 specific to pups. The elimination rate, accounting for all pup elimination routes, is: The removal rate constant for pups ( ,1 ) is defined analogously to females (Eq. (18)). The toxicokinetic equation (Eq. (41)) can be rewritten as the following differential equation for the PCB concentration 1, ( ) at time in a pup with a mother of age class during time period : d 1, d + ( ,1 + ,1 ) 1, = 1, 1, , + ,1 ,1 , = 2, … , 46 The body weight ratio between mother and pup is approximated as a mean value over the lactation period: Here, ̅ ( = 1, … , 46) are the mean body weights during the lactation period. The growth dilution rate ,1 is defined analogously to females (Eq. (20)) and is approximated as the mean growth dilution rate during a time period ( 0 ≤ ≤ d 1, d +̃, 1 1, = ̅ 1, , + ,1 ̅ ,1 , ̃1 , = 1, 1, , = 2, … , 46 , = , , The PCB concentration in the mother ( , ) during the lactation period ( = ) is obtained from Eq. (26) and inserted into Eq. (46), which yields the following first order ODE: The solution to Eq. (47) for a time period 0 ≤ ≤ ∆ within a year is: Equation (48) expresses the PCB concentration in a pup with a mother of age class at time during time period . The first term accounts for dietary uptake in the pup (proportional to dietary PCB concentration, characterized by constant ,1 ). The second term accounts for lactational transfer from mother to pup, reduced by elimination in the mother (exponential decrease, characterized by elimination rate ̃, ). The last term accounts for elimination in the pup (exponential decrease, characterized by elimination rate ̃, 1 ). Equation (48)is used to obtain the PCB concentration in pups at the end of different time periods , based on initial PCB concentrations , (0) at the start of the periods. The concentrations obtained at the end of a period is used as initial value for the next.
The initial PCB concentration in pups with mothers of age class + 1 (at the start of the lactation period year + 1) is assumed to equal the PCB concentration ̂0 , ( ) in fetuses with mothers of age class at the end of previous year ( ). The PCB concentration in pups with mothers of age class + 1 at the end of the lactation period ( = , = ∆ ) year + 1 is obtained from Eq. (48) as: The PCB concentration in pups with mothers of age class at the end of the delay period ( = , = ∆ ) year + 1 is obtained from Eq. (48) as: The PCB concentration in pups with mothers of age class at the end of the gestation period ( = , = ∆ ) year + 1 is obtained from Eq. (48) as: Equations (49), (50) and (51) are used to successively calculate PCB concentrations in pups at the end of the lactation period, the delay period and the gestation period for one year at a time.
Females that enter age class 2 are assumed to have an initial PCB concentration that equals the mean PCB concentration in one year old pups at the end of previous year (just before the occurrence of deaths and births). The mean PCB concentration in pups at the end of year + 1 is calculated as: Here, ( ) and ( ) is the survival rate and fertility of females in age class at the end of year , whereas ( − 1) is the number of females in age class at the end of year − 1 (after the occurrence of deaths). Hence, ( ) ( ) ( − 1) is the number of pups born by mothers of age class at the end of year , i.e. the number of pups with mothers of age class + 1 during year + 1.

Toxicodynamic Model
The toxicodynamic model describes how PCB body loads of fetuses, pups and females cause cumulation of damage, that results in hazard and reproductive stress and how these are translated into reduced survival and fertility (Fig. 6). The toxicodynamics are based on the threshold damage model (TDM), introduced by Ashauer et al. (2007). TDM is in turn based on DEBtox theory but adds cumulative effects and capacity to recover.

Fig. 6
The toxicodynamic model and its connection to other sub models. Text boxes represent sub modules. Free texts are inputs and outputs, with arrows showing relations to sub modules. The flowchart illustrates computations performed for year + 1, based on inputs from the previous year .
DEBtox assumes several physiological modes of action (pMoA) for toxicants in terms of which DEB parameters the toxicant stress affects (e.g. effects on growth, maturation and reproduction). Two PMoAs act directly on reproduction; reproductions costs (costs for maturation of females) and embryonic hazard (increased costs for somatic maintenance of embryos). Reproduction costs represent reduced efficiency in the utilization of energy for maturation of reproductive organs and offspring production, whereas embryonic hazard causes increased offspring mortality (an embryo dies if its somatic maintenance costs cannot be paid).
The current TD model follows the general principles of DEBtox in assuming that the toxic effects of PCBs on grey seals are considered through the lens of particular pMoAs affecting fundamental physiological processes like somatic maintenance. While our model is not based on DEB theory in the description of energetics and lifehistory traits, we interpret internal PCB concentration effects using a hazard model to account for mortality of adults, pups and fetuses, and damage to reproductive organs. In this context, we assume that decreased survival results from increased somatic maintenance costs (e.g., repair of lesions) due to cumulative damage causing hazard in pups and females. Decreased reproduction results from increased fetal mortality (fetal hazard) and reduced offspring production (reproductive stress). Instead of adopting a costs model, where negative reproductive effects of damaged reproductive organs result from increased costs for offspring production, we adopted a hazard model and considered reduced offspring production as a result of increased costs for maintenance of the reproductive apparatus due to cumulated damage (representing reproductive organs lesions etc.). In the hazard model, no offspring is produced if the maintenance costs of the reproductive apparatus cannot be paid. The costs model (Billoir, et al., 2007) was also explored to account for damaged reproductive organs, but found to be less consistent with empirical findings and thus not presented herein. Indirect effects on reproduction, such as decreased feeding rate and increased growth costs were neglected. In total, three different kinds of damage were considered: 1) damage to pups and females (causing hazard that decreases survival); 2) damage to reproductive apparatus (causing stress that reduces fertility); and 3) damage to fetuses (causing hazard that kills them and reduce the fertility of their mothers). Toxic effects and corresponding PMoAs are illustrated in Fig. 7. The hazard rate, describing the instantaneous probability to die, is defined as (Jager et al. 2011): Here, ℎ( ) is the hazard (the probability to die at time ) and ( ) is probability to survive until time . Hazard rate multiplied by a small time increment gives increase of hazard (probability to die during a time interval). The survival probability (the probability of the organism to survive until time ) is obtained from Eq. (53) as ( ) = 0 −ℎ( ) . The background survival (survival probability of animals with no exposure to toxicants) is 0 = − 0 , where 0 is the background mortality rate (mortality rate ( ) in absence of toxicant exposure). The probability of surviving from one age class ( − 1) to the next ( ) is given as: The hazard model (Billoir et al. 2007) accounts for fetal hazard by interpreting effects on reproduction rate as fetal mortality during gestation. Reproductive stress ( ), accounting for damaged reproductive apparatus at exposure time , is assumed to affect reproduction rate ( ) in a similar way (applying the hazard model): The fertility of age class is the number of offspring that reach age class 1 per female of age class :

Cumulation of Damage, Hazard and Stress
In order to calculate hazard and stress, a TDM similar to that established by Ashauer et al. (2007) was applied. The cumulation of all three types of damage (damage to pups/females, fetuses and reproductive organs) are described by the same governing equations, which only differ in values of model parameters. According to TDM, the damage rate increases linearly with the internal PCB concentration, whereas the recovery rate is proportional the current damage: Here, , ( ) is the damage at time (during time period ), , ( ) is the internal PCB concentration (mg/kg) at time , is the killing or stress rate constant (kg/mg·yr -1 ) of age class and is the recovery rate constant (yr -1 ) of age class . Notice that = 0, 1, represents fetuses, pups and females, respectively, whereas = 2, … , 46 represents different age classes of females. Hazard (or stress) occurs if the damage exceeds a critical threshold level: Here, ℎ , ( ) is the hazard at time (during time period ) and , is the damage threshold level. The operator [ ] + = max(0, ) is applied to ensure positive hazards. Equation (58) is a simplified way of expressing hazard compared to standard TDM (Ashauer et al. 2007) which applies a differential equation. In our model, both damage and hazard are nondimensional variables. Equations (57) and (58) are expressed in notations used for damage that causes decreased survival of females, pups and fetuses, but analogous equations are applied to describe reproductive damage ̃, and reproductive stress , , only differing in model parameters (̃, ̃ and ̃, ). In DEB theory, embryos and juveniles invest no energy in reproduction. In line with this, fetuses and pups are assumed to be unaffected by stress to the reproductive apparatus (̃0 = 1 = 0), which yields ̃0 , =̃1 , = 0 and 0, = 1, = 0. Damage to reproductive organs starts to cumulate when seals enter age class 2. The assumption is probably realistic since PCBs are hormone-like substances and should mainly affect females, with higher oestrogen levels than fetuses and pups.
Insertion of Eq. (59) into the toxicodynamic equation for damage rate (Eq. (57)) yields the following first order ODE: The solution to Eq. (61) is: Equation (62) can be used to calculate the cumulated damage , ( ) in fetuses ( = 0), pups ( = 1) or females ( = ) at a time during a time period within a year when the initial damage , (0) at the start of the period is known. We used Eq. (62) to successively calculate cumulated damage at the end of each time period ( = , , ) one year at a time. These included damages in pups and females (̂1 , ( ) and ̂( )), fetal damages ̂0 , ( ) and female reproductive damages ̃( ).
Obtained damage levels where used as initial values in the next calculation step, similar to how PCB concentrations were successively calculated. The damage level in females of age class + 1 at the start of year + 1 thus equals the damage level in females of age class at the end of previous year . The damage level in new-born pups with mothers of age class + 1 equals the damage level in fetuses with mothers of age class at the end of previous year. The initial damage in fetuses is zero. The initial damage in females of age class 2 at the start of a year was put to the mean damage in pups at the end of previous year, calculated similar to how concentrations were treated in Eq. (52): Here, ̂1 , +1 ( ) is the cumulated damage level in pups with mothers of age class + 1 at the end of the gestation period ( = ) year (the end of year t). Notice that calculations for a certain time period require knowledge of damage levels as well as PCB concentrations at the end of previous period.
Based on obtained damage levels, Eq. (58) was used to calculate hazards in pups and females ĥ ( ) at the end of each year , fetal hazards ĥ 0, ( ) at the end of each year and female reproductive stresses ̂( ) at the end of the delay period each year .

Fertilities and Survival Rates
Fertilities are affected by PCB exposure in two different ways; reduction due to reproductive stress and reduction due to fetal hazard. It was assumed that the state of females at the end of the delay period (just before gestation starts) determines the fertility reduction due to reproductive stress and that the state of fetuses at the end of the gestation period determines the fertility reduction due to fetal hazard.
The fertility of females in age class at the end of the delay period year was obtained from a hazard model (Eqs. (55) and (56)), accounting for decreased pregnancy rates due to reproductive organ lesions: The fertility of females in age class at the end of year was obtained by multiplying the fertility at the end of the delay period with a reduction factor accounting for fetal hazard ĥ 0, : Survival rates of pups and females are decreased by hazards ĥ . The survival rate for seals of age class at the end of year was thus obtained from Eq. (54) as: The ideal fertilities and survival rates ( 0 and 0 ) are vital rates for a population with maximal possible growth rate (in absence of toxicant exposure). Like Harding, et al. (2007), a full age-structured Leslie matrix model for Baltic grey seals was adopted, using annual time steps and 46 age classes. Since an insufficient number of males is unlikely to restrict growth of a seal population (Harding et al. 2007), only the size of the female population was considered. Individuals of age class = 1 (with an age of 0-1 year) are referred to as pups, whereas individuals of age class = 2, … , 46 are referred to as females. Unlike Harding, et al. (2007), the Leslie matrix elements are not constants. Fertilities ( ) and survival rates ( ) are linked to damage levels that have cumulated as a result of the history of internal PCB concentrations, which in turn depend on the history of PCB concentrations in the prey. Hence, the elements of the Leslie matrix are changing over time. It was assumed that all deaths occur at the end of a year and immediately after, all pups are born. Since only females that survive to the end of a year give birth to new pups, the first row of the Leslie matrix contains products of survival rates and fertilities, also a difference from Harding, et al. (2007). The Leslie matrix model is formulated as:

Leslie matrix model
The population vector ( ) includes the number of individuals in different age classes at the end of year . The Leslie matrix ( ) includes fertilities and survival rates at the end of year . Equation (68) is used to calculate the population vector at the end of a year, based on the population vector at the end of previous year. Since the population size is accounted for just after breeding, the model may be characterized as a post-breeding model. Ideal fertilities ( 0 ) and survival rates ( 0 ), corresponding to a grey seal population with maximal possible growth rate (in absence of toxicant exposure), were obtained from Harding, et al. (2007).
PCBs concentrations, damage levels and hazard/stress levels in fetuses, pups and females at the end of the different time periods within a year are successively calculated. Fertilities and survival rates are updated based on these computations and a new population vector is calculated, including the number of individuals in different age classes at the end of a year. Outputs from calculations for one year are used as inputs to calculations for the next year. The procedure is repeated for one year at a time and may be continued for as many years as desired.

Density dependence
Reduced fertilities due to density dependence are calculated in accordance with Caswell (2001) as: Here, is the ith row of , where is the effect of an individual in age class j on the fertility of an individual in age class i.
Reduced survival rates due to density dependence are similarly calculated as: Here, is the ith row of , where is the effect of an individual in age class j on the survival of an individual in age class i. To include density dependence in the TKTD population model, fertilities and survival rates in the Leslie matrix (Eq. (69)) are substituted by reduced values according to Eqs. (70) and (71). Unless otherwise stated, all performed analyses included density dependence.

Model Variables and Parameters
The TKTD population model simulates the temporal dynamics of a number of state variables for Baltic grey seals over their full lifespan (Table 1). We used the model to run simulations for the time period 1966-2015, where data on PCB levels in prey species were available. Segerstedt (2019) compiled data from previous studies of sum-PCB concentrations in the main prey items (cod, sprat and herring) of Baltic grey seals from all regions of the Baltic Sea during 1966-2015. The data are provided as lipid weight concentrations clustered into time intervals of 5 years (Fig. 8). The prey data illustrate that PCB levels peaked in the 1970s for all three fish species. We used the time clustered prey PCB data to linearly interpolate over each five-year period to obtain yearly PCB concentrations in prey, used as input to the TK model.
Since PCB levels in Baltic fish were low in the early 1960s (Bignert et al. 1998), all PCB and damage levels were put to zero at the start of year 1961. Prey concentrations were assumed to increase linearly from zero (in 1961) to reported levels in 1966. The total population size year 1961 has been estimated to 17 639 seals, including males and females (Harding and Härkönen 1999). The initial population size in the simulation (accounting for females exclusively) was put to 0 = 17 639/2 ≈ 8820.

Baltic grey seal biology and ecology
The TK model of temporal PCB dynamics incorporates aspects of grey seal life history, feeding habits and growth, as these have important effects on tissue PCB concentrations. The parameters defining the biology and ecology of grey seals in the model are detailed in Tables 2-5.

Grey seal growth and life history parameters
Adult grey seal males have an average weight of 233 kg and mature females have an average weight of 155 kg (Härkönen 2016). Females mature at an age of 3-5 years whereas males mature at around 6 years, although not socially mature until an age of 8 years (Jefferson et al. 2008). The maximum life span for Baltic grey seals is about 35-45 years, though only a minority of the seals become older than 24 years (HELCOM 2018). Baltic grey seals breed and mate in February to March (Jüssi et al. 2008). The fertilized egg is implanted in the uterus after a delay period of about 100 days. The succeeding gestation period of fetal growth to parturition lasts for about 240 days. Mature females give birth to at most one pup a year (Bergman 2007). The new-born pup has a weight of about 12 kg (Härkönen 2016). The pup is nursed for 15-18 days and weaning is abrupt (Härkönen 2016). During the lactation period, pups daily consume about 4 kg of breast milk (Iverson 1993). The milk contains about 52 % fat. At weaning, the pup has increased its weight to 37-48 kg, whereas the mother's weight loss is massive (Härkönen 2016). Grey seals breed on pack ice formations if they are available and otherwise on land (Jüssi et al. 2008). During two weeks in May and June, hundreds of seals gather at common sites to moult, without feeding (Hansson et al. 2017). Body growth parameters were estimated based on body size data for Baltic grey seals (primarily obtained from Härkönen (2016)) combined with calculations according to section 2.1. Life history and body growth parameters are summarised in Table 2. Table 2. Life history and body growth parameters

Grey seal diet variables and parameters
The prey consumption factor , the normalisation constant of the allometric relation that describes prey consumption rate ̇ in terms of body weight (Eq. (15)), was estimated from empirical data. Hansson, et al. (2017), estimated the daily mean prey consumption of Baltic grey seals to 4.5-5.0 kg, which corresponds to an annual per capita consumption of 1750 kg. It was here assumed that the largest females (with a weight of max = 160 kg) have a yearly fish consumption of ̅ max fish = 1800 kg. The prey consumption rate ̇, for the largest females (with body weight ) was estimated as the empirically obtained maximum yearly fish consumption of large grey seal females ( ̅ max fish ) divided by the length of the total feeding period during a year (∆ + ∆ ). The prey consumption factor was thus calculated as: The mean body lipid indices ̅ , for different prey species (Table 3) were obtained based on their fat content according to Fiskbasen (1996). Table 3. Dietary parameters for seals of different age classes. Estimated fraction of total prey biomass for different species consumed by Baltic grey seals and corresponding prey preference indices. Estimations of total consumed prey biomass are based on gut analyses of Baltic grey seal males and females (Lundström et al. 2010) and females (Tverin et al. 2019). Φ : Fraction of total prey biomass. : Prey preference index. ̅ , : Mean body lipid for prey species (Fiskbasen 1996). The mean sum-PCB concentration in the diet of a seal in age class was calculated as: The mean PCB concentration ̅ , ( ) in prey during year was obtained from the interpolated data of sum-PCB concentrations in the three prey species, collected by Segerstedt (2019) (Fig. 8).

Conversion of PCB concentrations
Empirical data of PCB concentrations in fish and seals are usually presented on a lipid weight basis, whereas the TKTD population model expresses PCBs in total body concentrations. To make comparisons, it is necessary to convert total body concentrations to lipid concentrations. The body fat index for seals of age class during time period is expressed as: Here, ̅ is the mean body weight (kg) of age class during time period and ̅ lip, is the corresponding body lipid content (kg). As a simplification it was assumed that all PCB in the body of a seal is bounded in fat tissues. The PCB lipid concentration , , can then be expressed in terms of the total body PCB concentration , as: , , = , / ̅ Eq. (76) was used to convert total body PCB concentrations to PCB lipid concentrations. For grey seals, most samples are collected in the autumn, at the middle of the gestation period when seals of all age classes are quite fat (Kauhala et al. 2017). As a simplification, it was assumed that mean body fat indices during sampling are the same for all age classes and remain constant throughout the gestation period; ̅ =̂. Hence, the PCB lipid concentrations at the end of year were calculated as: Analogously, the total body PCB concentration ̅ , ( ) of prey during year was calculated from the PCB lipid concentration ̅ , lip ( ) as: Here, ̅ , is the mean body lipid index for prey (Table 3). It was assumed that the mean body lipid index for Baltic grey seals during the gestation period is ̂= 30 % for all age classes (Iverson 2002).

Grey seal population parameters
Ideal vital rates ( 0 and 0 ) for a population with maximal possible reproduction were used as background values, representing fertilities and survival rates in absence of PCB exposure. These were obtained from Harding et al (2007), who reviewed life history data for grey seals from different parts of the world, derived vital rates for Baltic grey seals and established a Leslie matrix model, generating the observed population growth rate of = 1.075. They also derived optimal vital rates and calculated a maximum possible population growth rate of = 1.10. Values of vital rates are presented in Table 4. Notice that females younger than six years ( = 1 − 5) do not reproduce and that the youngest mature females ( = 6) have a fertility half that of adult females ( = 7, … , 45). In the ideal population, adult females obtain almost one pup a year, half of them female pups. Table 4. Age-specific fertilities and survival rates for Baltic grey seals. Realistic values ( , ) are values adopted in a population assessment by Harding, et al. (2007). Ideal values ( 0 , 0 ) correspond to a population with maximal possible growth rate (Harding et al. 2007). The formers are used as background values in the TKTD population model. Since seals are long-lived animals with high adult survival and low fecundity, populations typically consist of a large proportion of adults and the growth rate is low. This makes seal populations relatively resistant to short-term changes in the environment and the population dynamics are primarily governed by competition, predation, disease and availability of prey (Svensson et al. 2010). For long-lived K-strategic species, such as seals, density dependence is an important factor for longterm population dynamics (Kauhala et al. 2014). Seal populations show an age-specific response to increased population density with a rapid decrease in juvenile survival, decreased fecundity and a relatively constant adult survival, changing the age composition towards higher proportions of adults (Svensson et al. 2010). When a population approaches carrying capacity, the mortality is increased foremost among pups and sub-adults of 1-3 years. Hence, the most variable vital parameter in a seal population is pup survival. Pregnancy rates are lowest in the youngest females (4-5 years old), has a maximum during 6-24 years and then start to decline. However, females of 40 years may reproduce (Kauhala et al. 2014).
The effects of population density on survival and fertility of grey seals is modelled using age-specific survival and fertility density effect factors (Table 5). The carrying capacity was put to 100 000, since the Baltic grey seal population likely approached that size in the early 1900s (Harding and Härkönen 1999). It was assumed that adult survival and fertility are only dependent on the number of adults, not on the number of pups ( ,1 = ,1 = 0 , = 2, … , 46). Here, , and , are effects of an individual in age class on survival and fertility of an individual in age class . Since pups are more sensitive to harsh conditions than older animals, it was assumed that pup survival is affected by population density more than other transition rates ( 1, > , , 1, > , , , = 2, … , 46). It was also assumed that 1, = 1,1 ( = 2, … , 46), 1,1 = 2 2,2 , , = , = 2,2 ( , = 2, … , 46). The adult density factor 2,2 was adjusted such that the total population (females and males) approaches carrying capacity after long time under ideal conditions (no adverse effects from contaminants), yielding 2,2 = 2 · 10 −6 . Table 5. Density-dependence parameters. : Survival density effect (effect of an individual in age class on the survival of an individual in age class ). : Fertility density effect (effect of an individual in age class on the fertility of an individual in age class ). Fertility density effect of subadults and adults on adults 2 x10 -6

Toxicokinetics
Model parameters used in the TK model to describe temporal dynamics of PCB concentrations in seals throughout their lifespan are presented in Table 6, including rate constants that regulate PCB accumulation and depuration in grey seals, covering vertical transfer processes (placental and lactational transfer), dietary uptake, metabolic transformation, fecal egestion, and growth dilution. The derived parameter values are described below.

Assimilation efficiencies
In a bioaccumulation model for arctic ringed seals, Hickie, et al. (2005) adopted an assimilation efficiency of 90 %, both for uptake through diet and lactation. We adopted the same values for Baltic grey seals: Here, are the diet assimilation efficiency and is the lactation assimilation efficiency.

Vertical transfer rate constants
To estimate lactational transfer rate constants ( , ) and placental transfer rate constants ( , ), published data on milk transfer in lactating Canadian grey seals (Lang et al. 2011) was combined with published data on PCB transfer through breast milk in Scottish grey seals (Berghe et al. 2012). Lang, et al. (2011) presented body weight, body composition (fractions of water, protein and fat), daily milk output and milk composition for females and their pups at early and late lactation. A distinction was made between primiparous females (first-time breeders) and multiparous females (experienced breeders). Primiparous females had lower daily milk production, shorter lactation periods and smaller pups at weaning, compared to multiparous females. Berghe, et al. (2012) presented data on PCB concentrations on a lipid-weight basis for female-pup pairs at early and late lactation, including PCB levels in maternal outer and inner blubber, maternal blood, breast milk and pup blood. PCB concentrations increased significantly in all analysed tissues and milk between early and late lactation. In order to estimate lactational transfer rate constants from mentioned publications, it was assumed that the total PCB body burden in females decreases exponentially during lactation: Here, ( ) is the total PCB body burden of females at time during the lactation period (0 ≤ ≤ ∆ ). Lactational transfer rate constants were estimated as: In order to estimate total PCB body burdens ( ) at early lactation ( = ) and late lactation ( = ) from published data, the following assumptions were made: 1) all fat in a female body is either outer blubber or inner blubber and both layers are reduced during lactation; 2) the total amount of PCB in the outer blubber layer is constant throughout lactation and PCB concentration change is exclusively an effect of reduced blubber mass; 3) the inner blubber layer reduces linearly with time throughout lactation; 4) all inner blubber is consumed at the end of lactation; 5) the total PCB body burden of a female is reduced by the same amount as is transferred to her pup through breast milk during lactation; 6) milk PCB level increases linearly with time throughout lactation; 7) lost female body weight that is not blubber contains the PCB required to obtain the total PCB body burden according to assumption 5. The blood volume of grey seals is about 213 ml/kg and mammals have a blood density of 1060 kg/m 3 (Castellini and Mellish 2016). We assumed that about 10 % of a seal body is constituted by bone and estimated a blood weight fraction of 22.6 %. Assumptions 1-4 were used in calculations of outer/inner blubber ratios and assumptions 5-7 in calculations of total PCB body burdens.
The obtained lactational transfer rate constant for primiparous females ,6 = 4.36/yr corresponds to a 19 % reduction of the total PCB body burden during lactation (18 days), whereas the obtained lactational transfer rate constant for multiparous females , = 5.37/yr corresponds to a 23 % reduction. As a comparison, Addison and Brodie (1977) proposed that female grey seals lose about 15 % of their total PCB body burden during the nursing period.
In order to estimate placental transfer rate constants from mentioned publications, it was assumed that that the total PCB body burden in females during gestation decreases exponentially under constant body weight: Here, is total body PCB concentration in females at time ∆ during the gestation period (0 ≤ ∆ ≤ ∆ ). The placental transfer rate constant was estimated as: In order to estimate PCB concentrations before gestation (∆ = 0) and after gestation (∆ = ∆ ) from published data, the following assumptions were made: 1) blood PCB levels increase linearly in females and pups during lactation; 2) pup body and blubber weights increase exponentially from birth to late lactation; 3) the female outer blubber layer has the same weight at parturition and early lactation; 4) the blood/blubber PCB concentration ratio is the same for females and pups at parturition; 5) all PCB in females and pups at parturition is located in blubber; 6) the reduction of PCB body burden in females during gestation equals the total PCB body burden in pups at birth (dietary uptake and removal during gestation are neglected).
The obtained placental transfer rate constant , = 0.07/yr ( = 6, … , 46) corresponds to a 4 % loss of the total PCB body burden in females during gestation (247 days), which can be compared to the 19-23 % loss during lactation (18 days). This agrees with the general view that vertical PCB transfer in marine mammals is dominated by lactation, with only a minor contribution from placental transfer (Subramanian et al. 1988).

Removal rate constants
The removal rate constants , = , + , account for metabolic transformation and fecal egestion. According to Klanjscek et al. (2007), metabolic transformation of PCB in marine mammals is negligible, but we assume that a small fraction may be removed through metabolic transformation. In their bioaccumulation model for artic ringed seals, Hickie, et al. (2005) estimated elimination rate constants for 20 different PCB congeners. Linear regression was used to fit model outputs with observed contaminant concentrations in the seals. The estimated elimination rate constants ranged from 0.03/yr (PCB 153) to 2.5/yr (PCB 18), with 0.17/yr for sum-PCB. We considered these as removal rate constants, accounting for both metabolic transformation and fecal egestion. It was assumed that removal rate is proportional to metabolic rate (rate of energy spent per unit time) and that the former scales allometrically within and between seal species according to Kleiber's law. Metabolic rate thus increases with body mass according to a power law (Kleiber 1962): The normalisation constant is assumed to be similar for similar species of mammals. Larger animals thus metabolize PCB at a higher total rate than smaller animals, but smaller animals metabolize PCB faster per unit of body weight. The specific PCB removal rate is assumed to be proportional to the specific metabolic rate; ̇ ⁄ ∝ ⁄ ∝ −1 4 ⁄ . It then follows that the removal rate constant ( ) scales with body weight ( ) as ∝ −1/4 .
Adult ringed seals have a body weight of 50-70 kg (Jefferson et al. 2008). By assuming a mean body weight of = 60 kg, removal rate constants for ringed seals ( , ) can be transformed to removal rate constants for grey seals ( , ) of different age classes : Here, ̅ is the mean the body weight for age class during time period (where = 0 represents fetuses). The removal rate constant for sum-PCB in arctic ringed seals , = 0.17/yr (Hickie et al. 2005) was adopted to calculate values of removal rate constants for Baltic grey seals of different age classes, ranging from 0.1/yr (pregnant females) to 0.5/yr (fetuses).

Toxicodynamics
Model parameters used in the TD model to describe temporal dynamics of PCB effects in seals throughout their lifespan are presented in Table 7. The TD parameters define the rate constants, threshold levels, and tolerance levels that regulate the impact of PCB concentrations on hazard and stress in grey seals of different age classes. The derived parameter values are described below.
The toxicodynamic model includes three sets of parameters; reproductive stress parameters, fetal hazard parameters and hazard parameters related to survival of pups and adults. Each parameter set includes three type of constants; stress/killing rate constants, damage threshold levels and recovery rate constants. With respect to the many empirical observations of reproductive organ lesions in Baltic grey seals during the 1970s and the 1980s, stress to the reproductive apparatus is probably the most crucial path in which PCB affects vital rates. The corresponding reproductive stress parameters were estimated by a calibration procedure where pregnancy rates predicted by the TKTD population model were compared to reported pregnancy rates (Roos et al. 2012). For the other (less crucial) toxicodynamic parameters, values were chosen based on the reproductive stress parameters and on DEBtox parameters from minks (Desforges et al. 2017).

Calibration of Reproductive Stress Parameters
There are a number of publications, reporting historical data on pregnancy rates and prevalence of reproductive organ lesions associated with sterility (uterine occlusions, uterine stenosis and uterine leiomyoma) in Baltic grey seal females (Bergman 1999, Bäcklin et al. 2003, Bergman 2007, Bredhult et al. 2008, Roos et al. 2012, Bäcklin and Moraeus 2013 We compared this curve to pregnancy rates predicted by the TKTD population model, using temporal data on PCB concentrations in prey (Segerstedt 2019) as model input. It was assumed that reduced pregnancy is an effect of reproductive stress, caused by damage of the reproductive apparatus. In order to obtain the mean pregnancy rate ( ) from the TKTD population model, the mean fertility for age class 6-46 at the end of the delay period were multiplied by 2 (to include offspring of both sexes): The fertility rate at the end of the year was not used in the calculation, since this fertility also account for reproductive failure due to fetal hazard in pregnant females. The mean fertility of females in age class 6-46 at the end of the delay period year + 1 is obtained as: Here, ( ) is the number of females in age class at the end of year . Pregnancy rates, derived from the TKTD population model, were compared to observed pregnancy rates by visual inspection of graphs (Fig. 9). Reproductive stress parameters (̃, ̃,̃, ) were adjusted to fit model and empirical curve (Table 7).

Estimation of Toxicodynamic Model Parameters Related to Survival
The stress rate constant with respect to fetal hazard 0 was translated from DEBtox parameters for embryonal survival of minks (Desforges et al. 2017), see below (section 5.3.3). The corresponding damage threshold level ,0 was chosen arbitrarily as ten times smaller than the damage threshold level for reproductive stress in females since younger age classes are known to be more sensitive. The recovery rate constant for fetal hazard 0 was then scaled from the reproductive stress constants; 0 = ( 0̃2−46 ) ⁄̃2 −46 .
The killing rate constant with respect to pup survival 1 was translated from DEBtox parameters for kit survival of minks (Desforges et al. 2017), see below (section 5.3.3). The same killing rate constant was chosen for adults; 2−46 = 1 . The damage threshold level ,2−46 was adjusted to obtain good agreement between model and historical data for population size. The corresponding damage threshold level for pups was set to half of this value; ,1 = 0.5 ,2−46 . The recovery rate constants with respect to pup and female survival were scaled from reproductive stress constants; = (̃2 −46 ⁄ )̃2 −46 .

Transformation of DEBtox parameters
The toxicodynamic model is a threshold damage model (TDM), including stress rate constants, recovery rate constants and damage threshold levels. It is not a DEBtox model and does not use stress tolerance and NEC (no effect concentration) as parameters. Here, a rough and simple way to transform DEBtox parameters for one mammal (mink) into TDM parameters for another mammal (grey seal) is suggested. The method can be used as a way to obtain preliminary values before a more accurate parametrization is performed. As a demonstration, DEBtox parameters for mink are transformed into TDM parameters for grey seal. The obtained killing rate constant was used in the TKTD population model Baltic grey seals. Desforges, et al. (2017) developed a DEBtox model for mink, parameterized from published data on growth and reproduction in captive minks fed with PCB-contaminated fish. Embryos were exposed to PCB through the placenta throughout the gestation period (∆ = 42 days). Physical modes of action (PMoAs) were identified by comparing dose-response-curves from the model and reported data. Tolerance parameters and no effect concentrations were adjusted to fit model outputs with data. The model accurately predicted PCB accumulation and negative effects on embryo survival, kit survival, growth and development. DEBtox parameters for mink according to Table 7 were estimated (Desforges et al. 2017). Now, the transformation of these parameters is considered.
First, DEBtox parameters are transformed into TDM parameters (for the same species). Assuming no recovery ( = 0), the damage rate according to TDM (Eq. (57)) is: Assuming no initial damage and constant internal concentration , the cumulated damage after exposure time ∆ is: The hazard ℎ is obtained from Eq. (58) The hazard according to DEBtox theory only depends on the current toxicant concentration and is independent of exposure time: Comparing Eqs. (90) and (91) yields expressions for TDM parameters expressed in DEBtox parameters. The killing rate constant is expressed in terms of the tolerance concentration T and the exposure time ∆ as: The damage threshold level T is expressed in terms of the no effect concentration NEC and the tolerance concentration T as: An upper limit for the recovery rate constant may be estimated by assuming zero damage rate: The maximum recovery rate constant is thus: Since TDM includes cumulative effects, not considered in a DEBtox model, translation of DEBtox parameters into TDM parameters, requires knowledge of experimental exposure time ∆ . Mink DEBtox parameters were converted to TDM parameters by applying Eqs. (93) and (95) with exposure time ∆ = 42 days and are presented in Table 7. The TDM parameters for mink were then converted to TDM parameters for grey seal. According to Baas and Kooijman (2015) there is a strong correlation between specific metabolic rate (equated with specific somatic maintenance rate) and sensitivity to a toxicant (defined as having low NEC). This makes it possible to scale NEC between different animals. From analyses with different combinations of toxicants and animals, Baas & Kooijman (2015) found a strong negative correlation between specific somatic maintenance ( ) and no effect concentration (NEC): Here, and are toxicant-specific constants. If and NEC are known for two animals, exposed to a specified toxicant, and can be calculated. If is known for a third species, NEC can be estimated also for this one. The specific somatic maintenance for American mink (Neovison vison) and grey seal (Halichoerus grypus) are found in the Add-My-Pet database (Kooijman et With known value on the PCB-specific parameter , DEBtox data for mink can be converted to DEBtox data for grey seal: It was here assumed that = 1 and that damage threshold levels can be converted in the same way as no effect concentrations: It now follows from Eq. (93) that tolerance parameters are the same for different mammals exposed to a specified toxicant ( = ). By assuming equal exposure times in Eq. (92), also killing rate constants are obtained the same: = TDM parameters for mink were converted to TDM parameters for grey seal by applying Eqs. (95), (99) and (100).

Steady state
If prey PCB concentrations are constant over time, a stable state will finally be reached where PCB, damage, hazard and stress levels are stationary. If density dependence is neglected, also fertilities and survival rates stabilize at constant values and the population grows (or declines) at constant rate. Under these conditions, biomagnification factors and stable population growth rate can be defined. The mean biomagnification factor for all age classes is: Here ̅ lip ( ) is the mean PCB lipid concentration in seals of all age classes and ̅ lip is the mean PCB lipid concentration in the diet of all seals. The stable population growth rate is calculated as the dominant eigenvalue of the Leslie matrix: Here, ( ) are the eigenvalues of the Leslie matrix. Different levels of constant PCB lipid concentration in prey (0-10 mg/kg) were used as inputs to the TKTD population model and 100 years were simulated to investigate temporal changes in mean PCB lipid concentrations and population size for Baltic grey seals.

Sensitivity analysis
A sensitivity analysis investigates how uncertainty in input parameters causes uncertainty in population growth and can be used to identify parameters that are critical for population viability (Lacy et al. 2018). Assuming steady state conditions without density dependence (when growth rate is constant and nonzero), sensitivity analyses were performed for three toxicokinetic parameters (lactational transfer rate constants , , placental transfer rate constants , and removal rate constants , ) as well as for all toxicodynamic parameters (killing/stress rate constants, damage threshold levels and recovery rate constants).
A constant PCB lipid concentration in prey, generating positive stable growth rate over time, was used as input and the investigated parameters were varied one at time, whereas other parameters were held constant. The mean prey lipid concentration was set to 2 mg/kg, corresponding to Baltic levels of the early nineties. The relative value of a parameter ( ) is defined as the perturbed value ( ) divided by the value adopted in the model ( ); = / . Stable population growth rate was plotted versus perturbed relative value of investigated model parameter.